diff --git a/src/solvers/contact_solver.cpp b/src/solvers/contact_solver.cpp
index 05c27ae..dd8892c 100644
--- a/src/solvers/contact_solver.cpp
+++ b/src/solvers/contact_solver.cpp
@@ -1,64 +1,67 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 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 "contact_solver.hh"
 #include "logger.hh"
 #include <iomanip>
 #include <iostream>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 ContactSolver::ContactSolver(Model& model, const GridBase<Real>& surface,
                              Real tolerance)
     : model(model), tolerance(tolerance),
       surface_stddev(std::sqrt(surface.var())) {
 
   auto discretization = model.getBoundaryDiscretization();
   auto n_points = model.getTraction().getNbPoints();
 
   if (n_points != surface.dataSize())
     TAMAAS_EXCEPTION("Model size and surface size do not match!");
 
   this->surface.wrap(surface);
   _gap = allocateGrid<true, Real>(model.getType(), discretization);
   model.registerField("gap", _gap);
 }
 
 /* -------------------------------------------------------------------------- */
-
-void ContactSolver::printState(UInt iter, Real cost_f, Real error) {
-  // UInt precision = std::cout.precision();
+void ContactSolver::printState(UInt iter, Real cost_f, Real error) const {
   if (iter % dump_frequency)
     return;
+  logIteration(iter, cost_f, error);
+}
+
+/* -------------------------------------------------------------------------- */
+void ContactSolver::logIteration(UInt iter, Real cost_f, Real error) const {
   Logger().get(LogLevel::info)
       << std::setw(5) << iter << " " << std::setw(15) << std::scientific
       << cost_f << " " << std::setw(15) << error << std::endl
       << std::fixed;
 }
 
 /* -------------------------------------------------------------------------- */
 void ContactSolver::addFunctionalTerm(
     std::shared_ptr<functional::Functional> func) {
   functional.addFunctionalTerm(std::move(func));
 }
 
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
diff --git a/src/solvers/contact_solver.hh b/src/solvers/contact_solver.hh
index 3cb7f63..db90235 100644
--- a/src/solvers/contact_solver.hh
+++ b/src/solvers/contact_solver.hh
@@ -1,85 +1,87 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 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 CONTACT_SOLVER_HH
 #define CONTACT_SOLVER_HH
 /* -------------------------------------------------------------------------- */
 #include "meta_functional.hh"
 #include "model.hh"
 #include "tamaas.hh"
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 class ContactSolver {
 public:
   /// Constructor
   ContactSolver(Model& model, const GridBase<Real>& surface, Real tolerance);
   /// Destructor
   virtual ~ContactSolver() = default;
 
 public:
   /// Print state of solve
-  virtual void printState(UInt iter, Real cost_f, Real error);
+  virtual void printState(UInt iter, Real cost_f, Real error) const;
+  /// Log iteration info
+  virtual void logIteration(UInt iter, Real cost_f, Real error) const;
   /// Get maximum number of iterations
   auto getMaxIterations() const { return max_iterations; }
   /// Set maximum number of iterations
   void setMaxIterations(UInt n) { max_iterations = n; }
   /// Get dump_frequency
   auto getDumpFrequency() const { return dump_frequency; }
   /// Set dump_frequency
   void setDumpFrequency(UInt n) { dump_frequency = n; }
   /// Add term to functional
   void addFunctionalTerm(std::shared_ptr<functional::Functional> func);
   /// Returns functional object
   const functional::Functional& getFunctional() const { return functional; }
 
   // Virtual functions
 public:
   /// Solve for a mean traction vector
   virtual Real solve(std::vector<Real> /*load*/) {
     TAMAAS_EXCEPTION("solve() not implemented");
   }
 
   /// Solve for normal pressure
   virtual Real solve(Real load) { return solve(std::vector<Real>{load}); }
 
   /// Accessor for tolerance
   TAMAAS_ACCESSOR(tolerance, Real, Tolerance);
 
 public:
   GridBase<Real>& getSurface() { return surface; }
   Model& getModel() { return model; }
 
 protected:
   Model& model;
   GridBase<Real> surface;
   std::shared_ptr<GridBase<Real>> _gap = nullptr;
   functional::MetaFunctional functional;
 
   Real tolerance;
   UInt max_iterations = 1000;
   Real surface_stddev;
   UInt dump_frequency = 100;
 };
 
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
 #endif  // CONTACT_SOLVER_HH
diff --git a/src/solvers/polonsky_keer_rey.cpp b/src/solvers/polonsky_keer_rey.cpp
index 8a70d82..ec215dc 100644
--- a/src/solvers/polonsky_keer_rey.cpp
+++ b/src/solvers/polonsky_keer_rey.cpp
@@ -1,299 +1,302 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 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 "logger.hh"
 #include "loop.hh"
 #include "model_type.hh"
 #include <boost/preprocessor/seq.hpp>
 #include <iomanip>
 /* -------------------------------------------------------------------------- */
 
 namespace 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(
         std::make_shared<functional::ElasticFunctionalGap>(*integral_op,
                                                            this->surface));
     break;
   }
   case pressure: {
     model.getBEEngine().registerNeumann();
     primal = pressure_view.get();
     dual = gap_view.get();
     this->functional.addFunctionalTerm(
         std::make_shared<functional::ElasticFunctionalPressure>(*integral_op,
                                                                 this->surface));
     break;
   }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real PolonskyKeerRey::solve(std::vector<Real> target_v) {
   // Update in case of a surface change between solves
   this->surface_stddev = std::sqrt(this->surface.var());
 
   // Printing column headers
   Logger().get(LogLevel::info) << std::setw(5) << "Iter"
                                << " " << std::setw(15) << "Cost_f"
                                << " " << std::setw(15) << "Error" << '\n'
                                << std::fixed;
 
   const Real target = target_v.back();
   // 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;
 
   // Algorithm variables
   Real Gold = 1, error = 0;
   UInt n = 0;
   bool delta = false;
   *search_direction = 0;
 
   do {
     // Computing functional gradient
     functional.computeGradF(*primal, *dual);
 
     const 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
     const Real G = computeSquaredNorm(*dual);
 
     // Updating search direction (conjugate gradient)
     updateSearchDirection(delta * G / Gold);
     Gold = G;
 
     // Computing critical step
     const Real tau = computeCriticalStep(target);
 
     // Update primal variable
     delta = updatePrimal(tau);
 
     // We should scale to constraint
     if (constraint_type == variable_type)
       enforceMeanValue(target);
 
     error = computeError();
     const Real cost_f = functional.computeF(*primal, *dual);
     printState(n, cost_f, error);
   } while (error > this->tolerance and n++ < this->max_iterations);
 
   // Final update of dual variable
   functional.computeGradF(*primal, *dual);
 
+  // Log final iteration
+  logIteration(n, functional.computeF(*primal, *dual), error);
+
   enforceAdmissibleState();
 
   return error;
 }
 
 /* -------------------------------------------------------------------------- */
 void PolonskyKeerRey::enforceAdmissibleState() {
   /// Make dual admissible
   const 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;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * 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
   const Real sum = Loop::reduce<operation::plus>(
       [] CUDA_LAMBDA(Real & p, const Real& f) { return (p > 0) ? f : 0; },
       *primal, field);
   // Number of unsaturated points
   const 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);
 
   const Real rbar = meanOnUnsaturated(*projected_search_direction);
   if (constraint_type == variable_type)
     *projected_search_direction -= rbar;
   else {
     *projected_search_direction += 2 * target + rbar;
   }
 
   const 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);
   const 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) {
   const 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;
         }
         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();
 
   /// Complementarity error
   const 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.");
 
   const Real norm = [this]() {
     if (variable_type == pressure)
       return std::abs(primal->sum() * this->surface_stddev);
     return std::abs(dual->sum() * this->surface_stddev);
   }();
 
   return std::abs(error) / (norm * primal->getGlobalNbPoints());
 }
 
 /* -------------------------------------------------------------------------- */
 
 /**
  * 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();
 }
 
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */