diff --git a/src/solvers/beck_teboulle.cpp b/src/solvers/beck_teboulle.cpp
index 0c63308..494b5e3 100644
--- a/src/solvers/beck_teboulle.cpp
+++ b/src/solvers/beck_teboulle.cpp
@@ -1,111 +1,111 @@
 /**
  * @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 "beck_teboulle.hh"
 #include <iomanip>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 BeckTeboulle::BeckTeboulle(Model& model, const GridBase<Real>& surface,
                            Real tolerance, Real mu)
     : Kato(model, surface, tolerance, mu) {}
 
 /* -------------------------------------------------------------------------- */
 
 Real BeckTeboulle::solve(GridBase<Real>& g0) {
-  if (g0.getNbPoints() != comp) {
+  if (g0.getNbPoints() != pressure->getNbComponents()) {
     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 = solveTmpl<model_type::surface_1d>(g0);
     break;
   case model_type::surface_2d:
     cost = solveTmpl<model_type::surface_2d>(g0);
     break;
   default:
     break;
   }
 
   return cost;
 }
 
 template <model_type type>
 Real BeckTeboulle::solveTmpl(GridBase<Real>& g0) {
   // Printing column headers
   std::cout << std::setw(5) << "Iter"
             << " " << std::setw(15) << "Cost_f"
             << " " << std::setw(15) << "Error" << '\n'
             << std::fixed;
 
   constexpr UInt comp = model_type_traits<type>::components;
   Real cost = 0;
   UInt n = 0;
   *pressure = 0;
   GridBase<Real> z(*pressure);
   GridBase<Real> pressure_old(*pressure);
   Real t0 = 1;
   Real t1 = 1;
   Real m = 0;
   Real lipschitz = engine.getNeumannNorm();
 
   do {
     engine.solveNeumann(*pressure, *gap);
     addUniform<comp>(*gap, g0);
     *gap -= surfaceComp;
 
     z = *gap;
     z /= -lipschitz;
     *pressure += z;
     enforcePressureCoulomb<type>();
     z = *pressure;
 
     t1 = 1 + std::sqrt(4 * t0 * t0 + 1) / 2;
     m = 1 + (t0 - 1) / t1;
     z -= pressure_old;
     z *= m;
     z += pressure_old;
 
     t0 = t1;
     pressure_old = *pressure;
 
     cost = computeCost();
     printState(n, cost, cost);
   } while (cost > this->tolerance && n++ < this->max_iterations);
 
   computeFinalGap<type>();
 
   return cost;
 }
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
\ No newline at end of file
diff --git a/src/solvers/condat.cpp b/src/solvers/condat.cpp
index 5833dcb..0b42188 100644
--- a/src/solvers/condat.cpp
+++ b/src/solvers/condat.cpp
@@ -1,135 +1,135 @@
 /**
  * @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 "condat.hh"
 #include <iomanip>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 Condat::Condat(Model& model, const GridBase<Real>& surface,
                            Real tolerance, Real mu)
     : Kato(model, surface, tolerance, mu) {}
 
 /* -------------------------------------------------------------------------- */
 
 Real Condat::solve(GridBase<Real>& p0, Real grad_step) {
-  if (p0.getNbPoints() != comp) {
+  if (p0.getNbPoints() != pressure->getNbComponents()) {
     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, grad_step);
     break;
   case model_type::surface_2d:
     cost = solveTmpl<model_type::surface_2d>(p0, grad_step);
     break;
   default:
     break;
   }
 
   return cost;
 }
 
 template <model_type type>
 Real Condat::solveTmpl(GridBase<Real>& p0, Real grad_step) {
   // Printing column headers
   std::cout << std::setw(5) << "Iter"
             << " " << std::setw(15) << "Cost_f"
             << " " << std::setw(15) << "Error" << '\n'
             << std::fixed;
 
   constexpr UInt comp = model_type_traits<type>::components;
   Real cost = 0;
   UInt n = 0;
   pressure->uniformSetComponents(p0);
   GridBase<Real> pressure_old(*pressure); pressure_old = 0;
   Grid<Real, 1> q({1}, comp); q = 0;
 
   Real lipschitz = engine.getNeumannNorm();
   Real sigma = 2 * grad_step / lipschitz;
 
   do {
     updateGap<type>(sigma, grad_step, q);
     *pressure -= *gap;
     enforcePressureCoulomb<type>();
     updateLagrange<type>(q, pressure_old, p0);
     pressure_old = *pressure;
 
     cost = computeCost();
     printState(n, cost, cost);
   } while ((cost > this->tolerance || cost == 0.0) && n++ < this->max_iterations);
 
   computeFinalGap<type>();
 
   return cost;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Condat::updateGap(Real sigma, Real grad_step, GridBase<Real>& q) {
   constexpr UInt comp = model_type_traits<type>::components;
   engine.solveNeumann(*pressure, *gap);
 
   VectorProxy<Real, comp> _q = q(0);
   Vector<Real, comp> qt = _q;
   qt *= (1 - grad_step);
 
   Loop::stridedLoop(
     [sigma, grad_step, qt] CUDA_LAMBDA(VectorProxy<Real, comp>&& g,
                                VectorProxy<Real, comp>&& h) {
       g -= h;
       g *= sigma;
       g += qt;
     },
     *gap, surfaceComp);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Condat::updateLagrange(GridBase<Real>& q, GridBase<Real>& p_old,
                        GridBase<Real>& p0) {
   constexpr UInt comp = model_type_traits<type>::components;
 
   // GridBase<Real> corr(p0);
   Grid<Real, 1> corr({1}, comp);
   p_old *= -0.5;
   p_old += *pressure;
   computeMean<comp>(p_old, corr);
   corr *= 2;
   corr -= p0;
   q += corr;
 }
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
\ No newline at end of file
diff --git a/src/solvers/kato.cpp b/src/solvers/kato.cpp
index 37bf85b..92c0f58 100644
--- a/src/solvers/kato.cpp
+++ b/src/solvers/kato.cpp
@@ -1,361 +1,360 @@
 /**
  * @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, UInt proj_iter) {
-  if (p0.getNbPoints() != comp) {
+  if (p0.getNbPoints() != pressure->getNbComponents()) {
     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, proj_iter);
     break;
   case model_type::surface_2d:
     cost = solveTmpl<model_type::surface_2d>(p0, proj_iter);
     break;
   default:
     break;
   }
 
   return cost;
 }
 
 template <model_type type>
 Real Kato::solveTmpl(GridBase<Real>& p0, UInt proj_iter) {
   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, proj_iter);
     cost = computeCost();
     printState(n, cost, cost);
   } while (cost > this->tolerance && n++ < this->max_iterations);
 
   computeFinalGap<type>();
 
   return cost;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real Kato::solveRelaxed(GridBase<Real>& g0) {
-  if (g0.getNbPoints() != comp) {
+  if (g0.getNbPoints() != pressure->getNbComponents()) {
     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);
 
   computeFinalGap<type>();
 
   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, 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 proj_iter) {
   for (UInt i = 0; i < proj_iter; 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;
   corr *= -1;
   addUniform<comp>(*pressure, corr);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 void Kato::enforcePressureCoulomb() {
   constexpr UInt comp = model_type_traits<type>::components;
 
   Loop::stridedLoop(
     [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();
 
       // Projection normale au cône de friction
       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<operation::plus>(
     [] CUDA_LAMBDA(VectorProxy<Real, comp>&& f) -> Vector<Real, comp> {
       return f;
     },
     field);
 
   mean = _mean;
   mean /= N;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <UInt comp>
 void Kato::addUniform(GridBase<Real>& field, GridBase<Real>& vec) {
   VectorProxy<Real, comp> _vec(vec(0));
   field += _vec;
 }
 
 /* -------------------------------------------------------------------------- */
 
 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<operation::max>(
     [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(VectorProxy<Real, comp>&& p,
                              VectorProxy<Real, comp>&& g,
                              Real& lambda_,
                              Real& eta_,
                              Real& p_N_,
                              Real& p_C_) {
       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;
 
       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() {
   constexpr UInt comp = model_type_traits<type>::components;
   engine.solveNeumann(*pressure, *gap);
   *gap -= surfaceComp;
 
   Real g_N_min = Loop::stridedReduce<operation::min>(
     [] CUDA_LAMBDA(VectorProxy<Real, comp>&& g) {
       return g(comp - 1);
     },
     *gap);
 
   Grid<Real, 1> g_shift({comp}, 1);
   g_shift = 0;
   g_shift(comp - 1) = -g_N_min;
 
   *gap += surfaceComp;
   addUniform<comp>(*gap, g_shift);
 
   model.getDisplacement() = *gap;
 }
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
diff --git a/src/solvers/kato.hh b/src/solvers/kato.hh
index c36dadd..8e094b8 100644
--- a/src/solvers/kato.hh
+++ b/src/solvers/kato.hh
@@ -1,103 +1,102 @@
 /**
  * @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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __KATO_HH__
 #define __KATO_HH__
 /* -------------------------------------------------------------------------- */
 #include "contact_solver.hh"
 #include "meta_functional.hh"
 #include "model_type.hh"
 #include "static_types.hh"
 #include "tamaas.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 
 class Kato : public ContactSolver {
 public:
   /// Constructor
   Kato(Model& model, const GridBase<Real>& surface, Real tolerance, Real mu);
 
 public:
   /// Solve
   Real solve(GridBase<Real>& p0, UInt proj_iter);
   /// Solve relaxed problem
   Real solveRelaxed(GridBase<Real>& g0);
   /// Compute cost function
   Real computeCost();
 
 private:
   /// Template for solve function
   template <model_type type>
   Real solveTmpl(GridBase<Real>& p0, UInt proj_iter);
   /// Template for solveRelaxed function
   template <model_type type>
   Real solveRelaxedTmpl(GridBase<Real>& g0);
 
 protected:
   /// Creates surfaceComp form surface
   template <model_type type>
   void initSurfaceWithComponents();
   /// Project pressure on friction cone
   template <model_type type>
   void enforcePressureConstraints(GridBase<Real>& p0, UInt proj_iter);
   /// Project on C
   template <model_type type>
   void enforcePressureMean(GridBase<Real>& p0);
   /// Project on D
   template <model_type type>
   void enforcePressureCoulomb();
   /// Comupte mean value of field
   template <UInt comp>
   void computeMean(GridBase<Real>& field, GridBase<Real>& mean);
   /// Add vector to each point of field
   template <UInt comp>
   void addUniform(GridBase<Real>& field, GridBase<Real>& vec);
 
   /// Compute mean gap
   template <model_type type>
   Real computeBeta();
   /// Compute grids of dual and primal variables
   template <model_type type>
   void computeValuesForCost(Real beta, GridBase<Real>& lambda,
     GridBase<Real>& eta, GridBase<Real>& p_N, GridBase<Real>& p_C);
 
   /// Compute total displacement
   template <model_type type>
   void computeFinalGap();
 
 protected:
   BEEngine& engine;
   GridBase<Real>* gap = nullptr;
   GridBase<Real>* pressure = nullptr;
   Grid<Real, 1> surfaceComp;
   Real mu = 0;
-  UInt comp = 0; // number of components
   UInt N = 0; // number of points
 };
 
 __END_TAMAAS__
 
 #endif  // __KATO_HH__