diff --git a/src/bem_gigipol.cpp b/src/bem_gigipol.cpp
index deece03..63a6180 100644
--- a/src/bem_gigipol.cpp
+++ b/src/bem_gigipol.cpp
@@ -1,391 +1,389 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 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 <vector>
 #include "surface.hh"
 #include "bem_gigipol.hh"
 #include <iostream>
 #include <sstream>
 #include <fstream>
 #include <iomanip>
 #include <sstream>
 #include <cmath>
 /* -------------------------------------------------------------------------- */
 #define TIMER
 #include "surface_timer.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 Real BemGigipol::computeEquilibrium(Real epsilon,
                                  Real mean_displacement) {
 
 
     this->computeSpectralInfluenceOverDisplacement();
     this->surface_t = 0.;
     this->surface_r = 0.;
     this->surface_spectral_r = 0.;
     this->search_direction = 0.;
     this->pold = 0.;
     this->surface_displacements = 0.;
     Real delta = 0.;
     Real Gold = 1.;
     Real f = 1e300;
     Real fPrevious = 1e300;
     UInt n = surface.size();
     UInt size = n*n;
 
   this->true_displacements =0;
   this->true_displacements = mean_displacement;
   this->computeGaps();
   this->optimizeToMeanDisplacement(mean_displacement);
 
   convergence_iterations.clear();
   nb_iterations = 0;max_iterations=500;
 
   while(f > epsilon && nb_iterations < max_iterations) {
       fPrevious = f;
 
       this->functional->computeGradFU();
       this->search_direction= this->functional->getGradF();
       Real gbar = this->computeMeanPressuresInNonContact();
 
       this->search_direction  -= gbar;
       Real G = this->computeG();
       this->updateT(G,Gold,delta);
       Real tau = this->computeTau();
 
       this->old_displacements=this->true_displacements;
       Gold = G;
    #pragma omp parallel for
   for (unsigned int i = 0; i < size; ++i) {
 
               this->true_displacements(i)-=tau*this->surface_t(i);
       }
 
       //Projection on admissible space
     this->computeGaps();
 
 #pragma omp parallel for
   for (unsigned int i = 0; i < size; ++i) {
           if (this->gap(i)<0){
               this->true_displacements(i)=this->surface(i);}
       }
     this->computeGaps();
     delta = this->updateDisplacements(tau, mean_displacement);
 
   //    std::cout << "delta vaut "<<  delta << std::endl;
 
       //Projection on admissible space
       this->computeGaps();
       this->enforceMeanDisplacement(mean_displacement);
 
       this->computeGaps();
     this->computePressures();
       f=this->computeStoppingCriterion();
 
       if (nb_iterations % dump_freq == 0){
         Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions);
 
         std::cout << "G vaut "<< G<< std::endl;
         std::cout << "f vaut "<<  f << std::endl;
         std::cout << std::scientific << std::setprecision(10)
                   << nb_iterations << " "
                   << f << " " << f-fPrevious << " " << A
                   <<  std::endl;
 
       }
 
       convergence_iterations.push_back(f);
       ++nb_iterations;
 
 
 
 
 
     }
 
     this->computePressures();
 
 
     return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::computeStoppingCriterion() {
 
 
 
 
      Real crit=0.;
 
 
      UInt n = surface.size();
      UInt size = n*n;
 
    #pragma omp parallel for reduction(+:crit)
      for (UInt i = 0; i < size; ++i) {
         crit += std::abs(this->gap(i)*(this->surface_tractions(i)));
    ;
      }
 
 
 
   return crit ;
 }
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::computeG(){
   STARTTIMER("computeG");
 
   unsigned int n = surface.size();
   unsigned int size = n*n;
   Real res = 0.;
 #pragma omp parallel for reduction(+: res)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) > 0.) {
 
     Real val = this->search_direction(i);
     res += val*val;}
   }
 
   STOPTIMER("computeG");
   return res;
 }
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 Real BemGigipol::computeMeanPressuresInNonContact(){
   STARTTIMER("computeMeanPressuresInNonContact");
 
   unsigned int n = surface.size();
   unsigned int size = n*n;
   Real res = 0.;
   UInt nb_contact = 0;
 #pragma omp parallel for reduction(+: nb_contact,res)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) == 0. ) continue;
     ++nb_contact;
     res += this->search_direction(i);
   }
 
   res /= nb_contact;
   STOPTIMER("computeMeanPressuresInNonContact");
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 void BemGigipol::optimizeToMeanDisplacement(Real imposed_mean) {
   Real target_value = imposed_mean - SurfaceStatistics::computeAverage(surface, 0);
   UInt n = surface.size();
   UInt size = n * n;
 
   // Initial guesses for upper and lower bound
   Real step_min = -10;
   Real step_max = 10;
 
   // Gaps for upper and lower bound
   Real gap_min = positiveGapAverage(step_min);
   Real gap_max = positiveGapAverage(step_max);
 
   UInt max_expansion = 8;
 
   // Expand bounds if necessary
   for (UInt i = 0 ; gap_max < target_value && i < max_expansion ; i++, step_max *= 10)
     gap_max = positiveGapAverage(step_max);
 
   for (UInt i = 0 ; gap_min > target_value && i < max_expansion ; i++, step_min *= 10)
     gap_min = positiveGapAverage(step_min);
 
   Real g = 0.;
   Real epsilon = 1e-12;
   Real step = 0;
 
   while (fabs(g - target_value) > epsilon) {
     step = (step_min + step_max) / 2.;
     g = positiveGapAverage(step);
     if (g > target_value) step_max = step;
     else if (g < target_value) step_min = step;
     else {
       step_max = step;
       step_min = step;
     }
   }
 
   step = (step_min + step_max) / 2.;
 
 #pragma omp parallel for
   for (UInt i = 0 ; i < size ; i++) {
     gap(i) += step;
     if (gap(i) < 0) gap(i) = 0;
     true_displacements(i) = gap(i) + surface(i);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::positiveGapAverage(Real shift) {
   UInt n = surface.size();
 
   Real res = 0;
 #pragma omp parallel for reduction(+: res)
   for (UInt i = 0 ; i < n * n ; i++) {
     Real shifted_gap = gap(i) + shift;
     res += shifted_gap * (shifted_gap > 0);
   }
 
   return res / (n * n);
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 
 void BemGigipol::updateT(Real G, Real Gold, Real delta){
   STARTTIMER("updateT");
 
   unsigned int n = surface.size();
   unsigned int size = n*n;
   Real factor = delta*G/Gold;
 #pragma omp parallel for
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) == 0.) this->surface_t(i) = 0.;
     else{
       this->surface_t(i) *= factor;
       this->surface_t(i) += this->search_direction(i);
     }
   }
 
   STOPTIMER("updateT");
 }
 
 
 
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::computeTau() {
   STARTTIMER("computeOptimalStep");
 
   const UInt n = surface.size();
   const UInt size = n*n;
 
   this->applyInverseInfluenceFunctions(surface_t, surface_r);
 
   Real rbar = 0;
   UInt nb_contact = 0;
 #pragma omp parallel for reduction(+: nb_contact,rbar)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) == 0) continue;
     ++nb_contact;
     rbar += surface_r(i);
   }
   rbar /= nb_contact;
   surface_r -= rbar;
 
   Real tau_sum1 = 0.,tau_sum2 = 0.;
 #pragma omp parallel for reduction(+: tau_sum1, tau_sum2)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->gap(i) == 0) continue;
     tau_sum1 += this->search_direction(i)*surface_t(i);
     tau_sum2 += surface_r(i)*surface_t(i);
   }
   Real tau = tau_sum1/tau_sum2;
 
   STOPTIMER("computeTau");
   return tau;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::updateDisplacements(Real tau, Real mean_displacement) {
   STARTTIMER("updateDisplacements");
 
 
 
      unsigned int n = surface.size();
      unsigned int size = n*n;
 
 
 
      //compute number of interpenetration without contact
      UInt nb_iol = 0;
   #pragma omp parallel for reduction(+: nb_iol)
      for (unsigned int i = 0; i < size; ++i) {
        if ( this->search_direction(i) <0. && this->gap(i) == 0.){
 
          this->true_displacements(i) -= tau* this->search_direction(i);
          ++nb_iol;
        }
      }
 
      Real delta = 0;
      if (nb_iol > 0) delta = 0.;
      else            delta = 1.;
 
 
      return delta;
 
 
 
 
   STOPTIMER("updateDisplacements");
 }
 
 /* -------------------------------------------------------------------------- */
 void BemGigipol::enforceMeanDisplacement(Real mean_displacement) {
   STARTTIMER("enforceMeanDisplacement");
 
   UInt n = surface.size();
   UInt size = n*n;
   Real moyenne_surface=SurfaceStatistics::computeAverage(this->surface, 0);
   Real average = SurfaceStatistics::computeAverage(this->true_displacements, 0);
   Real factor = (mean_displacement-moyenne_surface) / (average-moyenne_surface);
 
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
  this->true_displacements(i) = factor*(this->true_displacements(i)-this->surface(i)) +this->surface(i);
   }
 
   STOPTIMER("enforceMeanDisplacement");
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigipol::computePressures() {
-  UInt n = surface.size();
-  UInt size = n*n;
   this->computeTractionsFromDisplacements();
   this->functional->computeGradFU();
   const Surface<Real> & gradF = this->functional->getGradF();
   Real min = SurfaceStatistics::computeMinimum(gradF);
   this->surface_tractions -= min;
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__