diff --git a/src/solvers/kato_saturated.cpp b/src/solvers/kato_saturated.cpp
index 7d41b4d..c74a7b2 100644
--- a/src/solvers/kato_saturated.cpp
+++ b/src/solvers/kato_saturated.cpp
@@ -1,159 +1,173 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "kato_saturated.hh"
 #include "logger.hh"
 #include <iomanip>
 #include <limits>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 KatoSaturated::KatoSaturated(Model& model, const GridBase<Real>& surface,
-                           Real tolerance, Real pmax)
+                             Real tolerance, Real pmax)
     : PolonskyKeerRey(model, surface, tolerance, PolonskyKeerRey::pressure,
                       PolonskyKeerRey::pressure),
       pmax(pmax) {}
 
 /* -------------------------------------------------------------------------- */
 Real KatoSaturated::meanOnUnsaturated(const GridBase<Real>& /*field*/) const {
   return 0;
 }
 
 Real KatoSaturated::computeSquaredNorm(const GridBase<Real>& /*field*/) const {
   return 1.;
 }
 
 void KatoSaturated::updateSearchDirection(Real /*factor*/) {
   *this->search_direction = *this->dual;
 }
 
 Real KatoSaturated::computeCriticalStep(Real /*target*/) {
   // integral_op->apply(*search_direction, *projected_search_direction);
   // Real num = search_direction->dot(*search_direction);
   // Real denum = projected_search_direction->dot(*search_direction);
   // return 0.1 * num / denum;
   return 1;
 }
 
 bool KatoSaturated::updatePrimal(Real step) {
   UInt na_num = Loop::reduce<operation::plus>(
       [step] CUDA_LAMBDA(Real & p, const Real& /*q*/, const Real& t) -> UInt {
         p -= step * t;  // Updating primal
         return 0;
         // if (p < 0)
         //   p = 0;                // Truncating neg values
         // if (p == 0 && q < 0) {  // If non-admissible state
         //   p -= step * q;
         //   return 1;
         // } else
         //   return 0;
       },
       *primal, *dual, *search_direction);
 
   return na_num == 0;
 }
 
 /* -------------------------------------------------------------------------- */
 Real KatoSaturated::computeError() {
   // We shift the gap by the minimum on unsaturated area
   Real shift = Loop::reduce<operation::min>(
       [&] CUDA_LAMBDA(const auto& p, const auto& d) {
         return (p < pmax) ? d : std::numeric_limits<Real>::max();
       },
       *primal, *dual);
   *dual -= shift;
 
   Real norm = 1;
 
   // Ignore points that are saturated
   Real error = Loop::reduce<operation::plus>(
       [&] CUDA_LAMBDA(auto p, auto d) { return (p < pmax) ? p * d : 0; },
       *primal, *dual);
 
   if (std::isnan(error))
     TAMAAS_EXCEPTION("Encountered NaN in complementarity error: this may be "
                      "caused by a contact area of a single node.");
 
   if (variable_type == pressure)
     norm = std::abs(primal->sum() * this->surface_stddev);
   else
     norm = std::abs(dual->sum() * this->surface_stddev);
 
   norm *= primal->getNbPoints();
   return std::abs(error) / norm;
 }
 
 /* -------------------------------------------------------------------------- */
 void KatoSaturated::enforceMeanValue(Real mean) {
   // We want to cancel the difference between saturated alpha + t and the
   // applied pressure
   *primal -= primal->mean();
   auto f = [&](Real scale) {
     Real sum = Loop::reduce<operation::plus>(
         [&] CUDA_LAMBDA(Real t) {
           t += scale;
 
           if (t > pmax)
             return pmax;
           else if (t < 0)
             return 0.;
           else
             return t;
         },
         *this->primal);
     sum /= this->primal->getNbPoints();
     return sum - mean;
   };
 
   if (pmax < mean)
     TAMAAS_EXCEPTION("Cannot find equilibrium");
 
   // Secant Newton on f
 
   // Initial points
   Real x_n_2 = 0., x_n_1 = pmax - primal->min(), x_n = 0.;
   Real f_n_2 = 0., f_n_1 = 0., f_n = 0.;
 
   // Secant loop
   do {
     f_n_2 = f(x_n_2);
     f_n_1 = f(x_n_1);
     x_n = x_n_1 - f_n_1 * (x_n_1 - x_n_2) / (f_n_1 - f_n_2);
     f_n = f(x_n);
     x_n_2 = x_n_1;
     x_n_1 = x_n;
   } while (std::abs(f_n / mean) > 1e-14);
 
   // Pressure update
   Loop::loop(
       [&] CUDA_LAMBDA(auto& t) {
         t += x_n;
         if (t > pmax)
           t = pmax;
         else if (t < 0)
           t = 0.;
       },
       *this->primal);
 }
 
+/* -------------------------------------------------------------------------- */
+void KatoSaturated::enforceAdmissibleState() {
+  /// Make dual admissible
+  Real shift = Loop::reduce<operation::min>(
+      [&] CUDA_LAMBDA(const auto& p, const auto& d) {
+        return (p < pmax) ? d : std::numeric_limits<Real>::max();
+      },
+      *primal, *dual);
+  *dual -= shift;
+
+  *displacement_view = *dual;
+  *displacement_view += this->surface;
+}
+
 }  // namespace tamaas
diff --git a/src/solvers/kato_saturated.hh b/src/solvers/kato_saturated.hh
index 73be942..2ddf3bd 100644
--- a/src/solvers/kato_saturated.hh
+++ b/src/solvers/kato_saturated.hh
@@ -1,60 +1,62 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef KATO_SATURATED_HH
 #define KATO_SATURATED_HH
 /* -------------------------------------------------------------------------- */
 #include "polonsky_keer_rey.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /// Polonsky-Keer algorithm with saturation of the pressure
 class KatoSaturated : public PolonskyKeerRey {
 public:
   /// Constructor
   KatoSaturated(Model& model, const GridBase<Real>& surface, Real tolerance,
                Real pmax);
   ~KatoSaturated() override = default;
 
 public:
   /// Mean on unsaturated constraint zone
   Real meanOnUnsaturated(const GridBase<Real>& field) const override;
   /// Compute squared norm
   Real computeSquaredNorm(const GridBase<Real>& field) const override;
   /// Update search direction
   void updateSearchDirection(Real factor) override;
   /// Compute critical step
   Real computeCriticalStep(Real target = 0) override;
   /// Update primal and check non-admissible state
   bool updatePrimal(Real step) override;
   /// Compute error/stopping criterion
   Real computeError() override;
   /// Enfore mean value constraint
   void enforceMeanValue(Real mean) override;
+  /// Enforce contact constraints on final state
+  void enforceAdmissibleState() override;
 
 protected:
   Real pmax;  ///< saturation pressure
 };
 
 }  // namespace tamaas
 #endif  // KATO_SATURATED_HH
diff --git a/src/solvers/polonsky_keer_rey.cpp b/src/solvers/polonsky_keer_rey.cpp
index 0693a62..18b891a 100644
--- a/src/solvers/polonsky_keer_rey.cpp
+++ b/src/solvers/polonsky_keer_rey.cpp
@@ -1,297 +1,303 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #include "polonsky_keer_rey.hh"
 #include "elastic_functional.hh"
 #include "fft_plan_manager.hh"
 #include "logger.hh"
 #include "loop.hh"
 #include "model_type.hh"
 #include <boost/preprocessor/seq.hpp>
 #include <iomanip>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 PolonskyKeerRey::PolonskyKeerRey(Model& model, const GridBase<Real>& surface,
                                  Real tolerance, type variable_type,
                                  type constraint_type)
     : ContactSolver(model, surface, tolerance), variable_type(variable_type),
       constraint_type(constraint_type) {
 
 #define SET_VIEWS(_, __, type)                                                 \
   case type: {                                                                 \
     setViews<type>();                                                          \
     break;                                                                     \
   }
 
   switch (model.getType()) {
     BOOST_PP_SEQ_FOR_EACH(SET_VIEWS, ~, TAMAAS_MODEL_TYPES);
   }
 
 #undef SET_VIEWS
 
   search_direction = allocateGrid<true, Real>(
       operation_type, model.getBoundaryDiscretization());
 
   projected_search_direction = allocateGrid<true, Real>(
 
       operation_type, model.getBoundaryDiscretization());
 
   switch (variable_type) {
   case gap: {
     model.getBEEngine().registerDirichlet();
     primal = gap_view.get();
     dual = pressure_view.get();
     this->functional.addFunctionalTerm(
         new functional::ElasticFunctionalGap(*integral_op, this->surface),
         true);
     break;
   }
   case pressure: {
     model.getBEEngine().registerNeumann();
     primal = pressure_view.get();
     dual = gap_view.get();
     this->functional.addFunctionalTerm(
         new functional::ElasticFunctionalPressure(*integral_op, this->surface),
         true);
     break;
   }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real PolonskyKeerRey::solve(std::vector<Real> target_v) {
   const Real target = target_v.back();
   Real G = 0, Gold = 1, error = 0, error_norm = 1;
   UInt n = 0;
   bool delta = false;
   *search_direction = 0;
   //*dual = 0;
 
   // Printing column headers
   Logger().get(LogLevel::info) << std::setw(5) << "Iter"
                                << " " << std::setw(15) << "Cost_f"
                                << " " << std::setw(15) << "Error" << '\n'
                                << std::fixed;
 
   // Setting uniform value if constraint
   if (constraint_type == variable_type && std::abs(primal->sum()) <= 0)
     *primal = target;
   else if (constraint_type == variable_type)
     *primal *= target / primal->mean();
   else if (constraint_type != variable_type)
     *primal = this->surface_stddev;
 
   do {
     // Computing functional gradient
     functional.computeGradF(*primal, *dual);
 
     Real dbar = meanOnUnsaturated(*dual);
     // Enforcing dual constraint via gradient
     if (constraint_type != variable_type) {
       *dual += 2 * target + dbar;
     } else {
       // Centering dual on primal > 0
       *dual -= dbar;
     }
 
     // Computing gradient norm
     G = computeSquaredNorm(*dual);
 
     // Updating search direction (conjugate gradient)
     updateSearchDirection(delta * G / Gold);
     Gold = G;
 
     // Computing critical step
     Real tau = computeCriticalStep(target);
 
     // Update primal variable
     delta = updatePrimal(tau);
 
     // We should scale to constraint
     if (constraint_type == variable_type)
       enforceMeanValue(target);
 
     error = computeError() / error_norm;
     Real cost_f = functional.computeF(*primal, *dual);
     printState(n, cost_f, error);
   } while (error > this->tolerance && n++ < this->max_iterations);
 
   // Final update of dual variable
   functional.computeGradF(*primal, *dual);
 
+  enforceAdmissibleState();
+
+  return error;
+}
+
+/* -------------------------------------------------------------------------- */
+void PolonskyKeerRey::enforceAdmissibleState() {
+  /// Make dual admissible
   Real dual_min = dual->min();
   *dual -= dual_min;
 
   // Primal is pressure: need to make sure gap is positive
   if (variable_type == pressure) {
     *displacement_view = *dual;
     *displacement_view += this->surface;
   }
   // Primal is gap: need to make sure dual is positive (pressure + adhesion)
   else {
     *displacement_view = *primal;
     *displacement_view += this->surface;
     integral_op->apply(*displacement_view, *pressure_view);
     *pressure_view -= dual_min;
   }
-
-  return error;
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * Computes \f$ \frac{1}{\mathrm{card}(\{p > 0\})} \sum_{\{p > 0\}}{f_i} \f$
  */
 Real PolonskyKeerRey::meanOnUnsaturated(const GridBase<Real>& field) const {
   // Sum on unsaturated
   Real sum = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f : 0; },
       *primal, field);
   // Number of unsaturated points
   UInt n_unsat = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(Real & p) -> UInt { return (p > 0); }, *primal);
 
   return sum / n_unsat;
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * Computes \f$ \sum_{\{p > 0\}}{f_i^2} \f$
  */
 Real PolonskyKeerRey::computeSquaredNorm(const GridBase<Real>& field) const {
   return Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f * f : 0; },
       *primal, field);
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * Computes \f$ \tau = \frac{ \sum_{\{p > 0\}}{q_i't_i} }{ \sum_{\{p > 0\}}{r_i'
  * t_i} } \f$
  */
 Real PolonskyKeerRey::computeCriticalStep(Real target) {
   integral_op->apply(*search_direction, *projected_search_direction);
 
   Real rbar = meanOnUnsaturated(*projected_search_direction);
   if (constraint_type == variable_type)
     *projected_search_direction -= rbar;
   else {
     *projected_search_direction += 2 * target + rbar;
   }
 
   Real num = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(const Real& p, const Real& q, const Real& t) {
         return (p > 0) ? q * t : 0;
       },
       *primal, *dual, *search_direction);
   Real denum = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(const Real& p, const Real& r, const Real& t) {
         return (p > 0) ? r * t : 0;
       },
       *primal, *projected_search_direction, *search_direction);
 
   return num / denum;
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * Update steps:
  * 1. \f$\mathbf{p} = \mathbf{p} - \tau \mathbf{t} \f$
  * 2. Truncate all \f$p\f$ negative
  * 3. For all points in \f$I_\mathrm{na} = \{p = 0 \land q < 0 \}\f$ do \f$p_i =
  * p_i - \tau q_i\f$
  */
 bool PolonskyKeerRey::updatePrimal(Real step) {
   UInt na_num = Loop::reduce<operation::plus>(
       [step] CUDA_LAMBDA(Real & p, const Real& q, const Real& t) -> UInt {
         p -= step * t;  // Updating primal
         if (p < 0)
           p = 0;                // Truncating neg values
         if (p == 0 && q < 0) {  // If non-admissible state
           p -= step * q;
           return 1;
         } else
           return 0;
       },
       *primal, *dual, *search_direction);
 
   return na_num == 0;
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * Error is based on \f$ \sum{p_i q_i} \f$
  */
 Real PolonskyKeerRey::computeError() {
   /// Making sure dual respects constraint
   *dual -= dual->min();
   Real norm = 1;
   Real error = primal->dot(*dual);
 
   if (std::isnan(error))
     TAMAAS_EXCEPTION("Encountered NaN in complementarity error: this may be "
                      "caused by a contact area of a single node.");
 
   if (variable_type == pressure)
     norm = std::abs(primal->sum() * this->surface_stddev);
   else
     norm = std::abs(dual->sum() * this->surface_stddev);
 
   norm *= primal->getNbPoints();
   return std::abs(error) / norm;
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * Do \f$\mathbf{t} = \mathbf{q}' + \delta \frac{R}{R_\mathrm{old}}\mathbf{t}
  * \f$
  */
 void PolonskyKeerRey::updateSearchDirection(Real factor) {
   Loop::loop(
       [factor] CUDA_LAMBDA(Real & p, Real & q, Real & t) {
         t = (p > 0) ? q + factor * t : 0;
       },
       *primal, *dual, *search_direction);
 }
 
 /* -------------------------------------------------------------------------- */
 void PolonskyKeerRey::enforceMeanValue(Real mean) {
-    *primal *= mean / primal->mean();
+  *primal *= mean / primal->mean();
 }
 
 /* -------------------------------------------------------------------------- */
 void PolonskyKeerRey::addFunctionalTerm(functional::Functional* func) {
   functional.addFunctionalTerm(func, false);
 }
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
diff --git a/src/solvers/polonsky_keer_rey.hh b/src/solvers/polonsky_keer_rey.hh
index 46af569..fa9538c 100644
--- a/src/solvers/polonsky_keer_rey.hh
+++ b/src/solvers/polonsky_keer_rey.hh
@@ -1,128 +1,130 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __POLONSKY_KEER_REY_HH__
 #define __POLONSKY_KEER_REY_HH__
 /* -------------------------------------------------------------------------- */
 #include "contact_solver.hh"
 #include "grid_view.hh"
 #include "meta_functional.hh"
 #include "westergaard.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 
 class PolonskyKeerRey : public ContactSolver {
 public:
   /// Types of algorithm (primal/dual) or constraint
   enum type { gap, pressure };
 
 public:
   /// Constructor
   PolonskyKeerRey(Model& model, const GridBase<Real>& surface, Real tolerance,
                   type variable_type, type constraint_type);
   ~PolonskyKeerRey() override = default;
 
 public:
   using ContactSolver::solve;
   /// Solve
   Real solve(std::vector<Real> target) override;
   /// Mean on unsaturated constraint zone
   virtual Real meanOnUnsaturated(const GridBase<Real>& field) const;
   /// Compute squared norm
   virtual Real computeSquaredNorm(const GridBase<Real>& field) const;
   /// Update search direction
   virtual void updateSearchDirection(Real factor);
   /// Compute critical step
   virtual Real computeCriticalStep(Real target = 0);
   /// Update primal and check non-admissible state
   virtual bool updatePrimal(Real step);
   /// Compute error/stopping criterion
   virtual Real computeError();
+  /// Enforce contact constraints on final state
+  virtual void enforceAdmissibleState();
   /// Add term to functional
   void addFunctionalTerm(functional::Functional* func);
   /// Enfore mean value constraint
   virtual void enforceMeanValue(Real mean);
 
 private:
   /// Set correct views on normal traction and gap
   template <model_type type>
   void setViews();
 
 protected:
   type variable_type, constraint_type;
   model_type operation_type;
   GridBase<Real>* primal = nullptr;  ///< non-owning pointer for primal varaible
   GridBase<Real>* dual = nullptr;    ///< non-owning pointer for dual variable
 
   /// CG search direction
   std::unique_ptr<GridBase<Real>> search_direction = nullptr;
   /// Projected CG search direction
   std::unique_ptr<GridBase<Real>> projected_search_direction = nullptr;
   /// View on normal pressure, gap, displacement
   std::unique_ptr<GridBase<Real>> pressure_view, gap_view,
       displacement_view = nullptr;
 
   /// Integral operator for gradient computation
   IntegralOperator* integral_op = nullptr;
 };
 
 /* -------------------------------------------------------------------------- */
 template <model_type mtype>
 void PolonskyKeerRey::setViews() {
   constexpr UInt dim = model_type_traits<mtype>::dimension;
   constexpr UInt bdim = model_type_traits<mtype>::boundary_dimension;
   constexpr UInt comp = model_type_traits<mtype>::components;
 
   pressure_view = std::unique_ptr<GridBase<Real>>{
       new GridView<Grid, Real, bdim, bdim>(model.getTraction(), {}, comp - 1)};
   gap_view = std::unique_ptr<GridBase<Real>>{
       new GridView<Grid, Real, bdim, bdim>(*this->_gap, {}, comp - 1)};
   displacement_view =
       std::unique_ptr<GridBase<Real>>{new GridView<Grid, Real, dim, bdim>(
           model.getDisplacement(), model_type_traits<mtype>::indices,
           comp - 1)};
 
 #define WESTERGAARD(type, kind, desc)                                          \
   this->model.registerIntegralOperator<                                        \
       Westergaard<model_type::type, IntegralOperator::kind>>(desc)
 
   // I noz is fugly - TODO factory?
   if (bdim == 1) {
     operation_type = model_type::basic_1d;
     if (variable_type == gap)
       integral_op = WESTERGAARD(basic_1d, dirichlet, "westergaard_dirichlet");
     else
       integral_op = WESTERGAARD(basic_1d, neumann, "westergaard_neumann");
   } else if (bdim == 2) {
     operation_type = model_type::basic_2d;
     if (variable_type == gap)
       integral_op = WESTERGAARD(basic_2d, dirichlet, "westergaard_dirichlet");
     else
       integral_op = WESTERGAARD(basic_2d, neumann, "westergaard_neumann");
   }
 
 #undef WESTERGAARD
 }
 
 __END_TAMAAS__
 
 #endif  // __POLONSKY_KEER_REY_HH__
diff --git a/tests/test_saturated_pressure.py b/tests/test_saturated_pressure.py
index 0fa9edc..d3a4acb 100644
--- a/tests/test_saturated_pressure.py
+++ b/tests/test_saturated_pressure.py
@@ -1,67 +1,66 @@
 # -*- coding: utf-8 -*-
 # @file
 # @section LICENSE
 #
 # Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
 # Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
 #
 # This program is free software: you can redistribute it and/or modify
 # it under the terms of the GNU Affero General Public License as published
 # by the Free Software Foundation, either version 3 of the License, or
 # (at your option) any later version.
 #
 # This program 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 Affero General Public License for more details.
 #
 # You should have received a copy of the GNU Affero General Public License
 # along with this program.  If not, see <https://www.gnu.org/licenses/>.
 
 from __future__ import print_function, division
 
 import numpy as np
 import tamaas as tm
 
 
 def test_saturated_pressure(tamaas_fixture):
     grid_size = 256
     load = 0.06
     p_sat = 0.4
 
     iso = tm.Isopowerlaw2D()
     iso.q0 = 4
     iso.q1 = 4
     iso.q2 = 16
     iso.hurst = 0.8
 
     sg = tm.SurfaceGeneratorFilter2D()
     sg.random_seed = 2
     sg.setSizes([grid_size, grid_size])
     sg.setSpectrum(iso)
     surface = sg.buildSurface()
     surface *= 0.01 / iso.rmsHeights()
     surface -= np.max(surface)
 
     model = tm.ModelFactory.createModel(tm.model_type.basic_2d,
                                         [1., 1.], [grid_size, grid_size])
     model.E = 1.
     model.nu = 0
 
     solver = tm.KatoSaturated(model, surface, 1e-12, p_sat)
     solver.setMaxIterations(6000)
     assert solver.solve(load) < 1e-12
 
     tractions = model['traction']
 
     mean_pressure_error = abs(np.mean(tractions)-load)/load
     max_pressure_error = np.max(tractions)-p_sat
 
     gaps = model['gap']
     gaps[np.where(tractions == p_sat)] = 0.
-    # press = tractions
-    # orthogonality_error = np.max(press*gaps) / (press.max() * gaps.max())
+    press = tractions
+    orthogonality_error = np.max(press*gaps) / (press.max() * gaps.max())
 
-    # assert mean_pressure_error < 1e-12 and max_pressure_error < 1e-12 and \
-    #     orthogonality_error < 2e-10
-    assert mean_pressure_error < 1e-12 and max_pressure_error < 1e-12
+    assert mean_pressure_error < 1e-12 and max_pressure_error < 1e-12 and \
+        orthogonality_error < 4e-8