diff --git a/src/adhesion_functional.cpp b/src/adhesion_functional.cpp
index 723e1de..c2dbf6f 100644
--- a/src/adhesion_functional.cpp
+++ b/src/adhesion_functional.cpp
@@ -1,132 +1,141 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@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 "adhesion_functional.hh"
 #include <cmath>
 
 /* -------------------------------------------------------------------------- */
 
 AdhesionFunctional::AdhesionFunctional(BemFFTBase & bem):
   ElasticEnergyFunctional(bem),
   rho(1.),
   surface_energy(1.),
   exponential_term(bem.getSurface().size(), bem.getSurface().getL()),
   spectral_gap(bem.getSurface().size(), bem.getSurface().getL())
 {}
 
 AdhesionFunctional::~AdhesionFunctional()
 {}
 
 void AdhesionFunctional::computeF() {
 
 
     Real res = 0;
 
 
     // this->computeExponentialTerm();
 
-    Surface<Real> & gap = const_cast<Surface<Real> &>(bem.getTrueDisplacements());
+    Surface<Real> & gap = const_cast<Surface<Real> &>(bem.getGap());
     gap-=SurfaceStatistics::computeAverage(  gap, 0);
 
 
 
     const Surface<Real> & influence_coefficients = this->bem.getSpectralInfluenceOverDisplacement();
 
     const UInt n = this->exponential_term.size();
     const UInt size = n * n;
 
     /// Computation of elastic energy gradient
     gap.FFTTransform(this->spectral_gap, bem.getNumberOfThreads());
 
     this->spectral_gap(0) = 0;
   #pragma omp parallel for
     for (UInt i = 1 ; i < size ; ++i){
 
       this->spectral_gap(i) /= influence_coefficients(i);
     }
+
+    // this->spectral_gap(n/2,n/2) = 0.;//1./d;
+    //i == 0 j == n/2
+
+    // this->spectral_gap(0,n/2) = 0.;//1./d;
+    //i == n/2 j == 0
+   //  this->spectral_gap(n/2,0) = 0.;//1./d;
+
+
     this->spectral_gap.FFTITransform(this->gradF, bem.getNumberOfThreads());
 
 
     for (UInt i = 0; i < size; ++i) {
 
       res += 0.5*this->gradF(i)*gap(i);
 
       // res +=
       //   this->surface_tractions[i].real()
       //   *(surface_displacements[i].real() - surface[i].real());
     }
 
     this->F = res;
 }
 
 void AdhesionFunctional::computeGradF() {
   // this->computeExponentialTerm();
 
-  Surface<Real> & dep = const_cast<Surface<Real> &>(bem.getTrueDisplacements());
+  Surface<Real> & gap =  const_cast<Surface<Real> &>(bem.getGap());
   const Surface<Real> & influence_coefficients = this->bem.getSpectralInfluenceOverDisplacement();
 
   const UInt n = this->exponential_term.size();
   const UInt size = n * n;
 
   /// Computation of elastic energy gradient
-  dep.FFTTransform(this->spectral_gap, bem.getNumberOfThreads());
+  gap.FFTTransform(this->spectral_gap, bem.getNumberOfThreads());
   this->spectral_gap(0) = 0;
 #pragma omp parallel for
   for (UInt i = 1 ; i < size ; ++i){
 
     this->spectral_gap(i) /= influence_coefficients(i);
 
 
 
 }
   this->spectral_gap.FFTITransform(this->gradF, bem.getNumberOfThreads());
 
   /// Computation of adhesion energy gradient
 
   this->computeExponentialTerm();
 
  #pragma omp parallel for
    for (UInt i = 0 ; i < size ; ++i){
 
     this->gradF(i) += this->exponential_term(i) / this->rho ;
 
 
    }
 }
 
 void AdhesionFunctional::computeExponentialTerm() {
   const Surface<Real> & gap = bem.getGap();
   const UInt n = gap.size();
   const UInt size = n * n;
 
 #pragma omp parallel for
   for (UInt i = 0 ; i < size ; ++i) {
     this->exponential_term(i) = this->surface_energy * exp(- gap(i) / this->rho);
 
 
   }
 }
diff --git a/src/bem_gigi.cpp b/src/bem_gigi.cpp
index a932e26..22ac84d 100644
--- a/src/bem_gigi.cpp
+++ b/src/bem_gigi.cpp
@@ -1,459 +1,463 @@
 /**
  *
  * @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_gigi.hh"
 #include <iostream>
 #include <sstream>
 #include <fstream>
 #include <iomanip>
 #include <sstream>
 #include <cmath>
 /* -------------------------------------------------------------------------- */
 #define TIMER
 #include "surface_timer.hh"
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeEquilibrium(Real epsilon,
                                  Real mean_displacement) {
 
   this->computeSpectralInfluenceOverDisplacement();
 
   Real Rold = 1.;
   Real Gold=1.;
   Real f = 1e300;
 
   this->search_direction = 0.;
   this->true_displacements =0;
   this->true_displacements = mean_displacement;
    this->computeGaps();
    this->optimizeToMeanDisplacement(mean_displacement);
    this->computeGaps();
 
 
 
 
   std::cout << "moyenne deplacement "<< SurfaceStatistics::computeAverage(this->true_displacements, 0)<< std::endl;
 
 
   convergence_iterations.clear();
   nb_iterations = 0;
- UInt crit2=50;
+ UInt crit2=100;
 
 
-   max_iterations =50;
+   max_iterations =100;
 
 
   while (f> epsilon && nb_iterations++ < max_iterations) {
 
     this->computeGaps();
     this->functional->computeGradF();
       Surface<Real> & gradF = this->functional->getGradF();
 
 
     this->functional->computeF();
     Real F=this->functional->getF();
      std::cout << "fonctionnal value "<< F << std::endl;
 
     Real R = this->functional->computeGradientNorm();
+
+    std::cout << "norme gradient "<< R << std::endl;
+
+
     this->updateSearchDirection(R / Rold);
     Rold = R;
     Real alpha = this->computeOptimalStep();
     this->old_displacements=this->true_displacements;
     this->updateDisplacements(alpha);
 
     this->computeGaps();
     this->optimizeToMeanDisplacement(mean_displacement); //espace admissible
+    this->computeGaps();
 
-
-    this->computePressures(mean_displacement); //p=K-1*u
+   // this->computePressures(mean_displacement); //p=K-1*u
 
     Real G=this->computeG();
 
-
+G=0.;
 
      // std::cout << "G vaut "<< G<< std::endl;
    UInt inner_loop_cpt=0;
 
        //Inner loop
-    while(G > 0 && inner_loop_cpt < crit2 ) {
+    while(G > 0. && inner_loop_cpt < crit2 ) {
       ++inner_loop_cpt;
-
+       this->computeGaps();
       this->functional->computeGradF();
       this->updateSearchDirection(R / Rold);
 
       Real tau = this->computeTau();
       Gold = G;
       this->updateDisplacementsB(tau);
       this->computeGaps();
       this->optimizeToMeanDisplacement(mean_displacement); //espace admissible
       this->computePressures(mean_displacement); //p=K-1*u
+      this->computeGaps();
       G = this->computeG();
            std::cout << "G vaut "<< G<< std::endl;
 
       }
 
     f = computeStoppingCriterion();
     if (nb_iterations % dump_freq == 0) {
       std::cout << std::scientific << std::setprecision(10)
                 << nb_iterations << " "
                 << f << std::fixed << std::endl;
     }
 
     convergence_iterations.push_back(f);
     // if(nb_iterations == 2) return 0;
   }
 
 
 
-    this->computePressures(mean_displacement);
- // this->computeTruePressures(mean_displacement);
+  //  this->computePressures(mean_displacement);
 
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeStoppingCriterion() {
 
 
 
 
      Real crit=0.;
      Real disp_norm = 0.;
 
      UInt n = surface.size();
      UInt size = n*n;
 
    #pragma omp parallel for reduction(+:crit, disp_norm)
      for (UInt i = 0; i < size; ++i) {
        crit += ( this->old_displacements(i) - true_displacements(i))*(
    this->old_displacements(i) - true_displacements(i));
        disp_norm += (true_displacements(i)*true_displacements(i));
      }
 
 
 
   return crit / disp_norm;
 }
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::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->surface_tractions(i);
     res += val*val;}
   }
 
   STOPTIMER("computeG");
   return res;
 }
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::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 BemGigi::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 BemGigi::updateSearchDirection(Real factor) {
   STARTTIMER("updateSearchDirection");
 
   UInt n = surface.size();
   UInt size = n*n;
   Surface<Real> & gradF = this->functional->getGradF();
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
-    //this->search_direction(i) *= factor;
+
     this->search_direction(i) = gradF(i);
   }
 
   STOPTIMER("updateSearchDirection");
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeOptimalStep() {
   STARTTIMER("computeOptimalStep");
 
   search_direction.FFTTransform(spectral_search_direction, nthreads);
 
   UInt n = surface.size();
   UInt size = n*n;
   spectral_search_direction(0) = 0;
 #pragma omp parallel for
   for (UInt i = 1; i < size; ++i) {
     spectral_search_direction(i) /= this->surface_spectral_influence_disp(i);
   }
 
   // /!\ does not contain the spectral search direction anymore
   spectral_search_direction.FFTITransform(nthreads);
  // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
  // spectral_search_direction -= average;
 
   Surface<Real> & gradF = this->functional->getGradF();
 
   Real numerator = 0., denominator = 0.;
 
 #pragma omp parallel for reduction(+: numerator, denominator)
   for (UInt i = 0; i < size; ++i) {
     numerator += gradF(i) * search_direction(i);
     denominator += spectral_search_direction(i).real() * search_direction(i);
   }
 
   Real alpha = numerator / denominator;
 
   STOPTIMER("computeOptimalStep");
   return alpha;
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeTau() {
   STARTTIMER("computeOptimalStep");
 
   search_direction.FFTTransform(spectral_search_direction, nthreads);
 
   UInt n = surface.size();
   UInt size = n*n;
   spectral_search_direction(0) = 0;
 #pragma omp parallel for
   for (UInt i = 1; i < size; ++i) {
     spectral_search_direction(i) /= this->surface_spectral_influence_disp(i);
   }
 
   // /!\ does not contain the spectral search direction anymore
   spectral_search_direction.FFTITransform(nthreads);
  // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
  // spectral_search_direction -= average;
 
   Surface<Real> & gradF = this->functional->getGradF();
 
   Real numerator = 0., denominator = 0.;
 
 #pragma omp parallel for reduction(+: numerator, denominator)
   for (UInt i = 0; i < size; ++i) {
 if (this->gap(i) > 0) {
     numerator += gradF(i) * search_direction(i);
     denominator += spectral_search_direction(i).real() * search_direction(i);}
   }
 
   Real tau = numerator / denominator;
 
   STOPTIMER("computeOptimalStep");
   return tau;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::updateDisplacements(Real alpha) {
   STARTTIMER("updateDisplacements");
 
   UInt n = surface.size();
   UInt size = n*n;
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     this->true_displacements(i) -= alpha*this->search_direction(i);
   }
 
   STOPTIMER("updateDisplacements");
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::updateDisplacementsB(Real alpha) {
   STARTTIMER("updateDisplacements");
 
   UInt n = surface.size();
   UInt size = n*n;
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
       if (this->gap(i)>0){
      this->true_displacements(i) -= alpha*this->search_direction(i);
 
       }
   }
 
   STOPTIMER("updateDisplacements");
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::emptyOverlap() {
   STARTTIMER("emptyoverlap");
 
   UInt n = surface.size();
   UInt size = n*n;
 
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     if (gap(i) < 0)
       this->true_displacements(i) = this->surface(i);
   }
 
   STOPTIMER("emptyoverlap");
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::enforceMeanDisplacement(Real mean_displacement) {
   STARTTIMER("enforceMeanDisplacement");
 
   UInt n = surface.size();
   UInt size = n*n;
   Real average = SurfaceStatistics::computeAverage(this->true_displacements, 0);
   Real factor = mean_displacement / average;
 
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     this->true_displacements(i) *= factor;
   }
 
   STOPTIMER("enforceMeanDisplacement");
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::computePressures(Real mean_displacement) {
   UInt n = surface.size();
   UInt size = n*n;
 
   Real moy_surface=SurfaceStatistics::computeAverage(this->surface, 0);
 
   this->true_displacements.FFTTransform(this->surface_spectral_tractions, nthreads);
   this->surface_spectral_tractions(0) = 0;
 
 #pragma omp parallel for
   for (UInt i = 1 ; i < size ; ++i) {
     this->surface_spectral_tractions(i) /= this->surface_spectral_influence_disp(i);
   }
 
   this->surface_spectral_tractions.FFTITransform(this->surface_tractions, nthreads);
   Real true_pressure = 0.;
   for (UInt i = 0 ; i < n * n ; i++) {
       true_pressure+=this->surface_tractions(i)*(this->true_displacements(i)-mean_displacement-this->surface(i));
 
   }
 
   true_pressure/=(moy_surface*size-mean_displacement*size);
-
+     std::cout << "true_pressure "<< true_pressure << std::endl;
    this->surface_tractions += true_pressure;
 
   //this->surface_tractions=  this->old_displacements - true_displacements;
 
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::computeTruePressures(Real mean_displacement) {
   this->computeGaps();
  this->functional->computeGradF();
     UInt n = gap.size();
 
   //Orthogonalite
   for (UInt i = 0 ; i < n * n ; i++) {
       if (this->gap(i) > 0){
 
              this->surface_tractions(i)=0;
 
 
       }
   }
 
 
 
 
 }