diff --git a/src/bem_gigi.cpp b/src/bem_gigi.cpp
index fdd6f89..2957f35 100644
--- a/src/bem_gigi.cpp
+++ b/src/bem_gigi.cpp
@@ -1,362 +1,321 @@
 /**
  *
  * @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"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 Real BemGigi::computeEquilibrium(Real epsilon,
                                  Real mean_displacement) {
 
   this->computeSpectralInfluenceOverDisplacement();
 
   Real Rold = 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();
 
 
 
   convergence_iterations.clear();
   nb_iterations = 0;
   max_iterations =50;
    Real R=0;
 
   while (f> epsilon && nb_iterations++ < max_iterations) {
 
     this->computeGaps();
     this->functional->computeGradFU();
     this->updateSearchDirection(R / Rold);
 
     Rold = R;
     Real alpha = this->computeOptimalStep();
         std::cout << "alpha vaut "<< alpha<< std::endl;
 
     this->old_displacements=this->true_displacements;
     this->updateDisplacements(alpha);
     this->computeGaps();
 
     this->optimizeToMeanDisplacement(mean_displacement); //espace admissible
        this->computeGaps();
   this->computePressures(mean_displacement);
     f=computeStoppingCriterion();
 
 }
   this->computePressures(mean_displacement);
 
 
 
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeStoppingCriterion() {
 
 
     Real rho = this->functional->getParameter("rho");
     Real surface_energy = this->functional->getParameter("surface_energy");
 
 
      Real crit=0.;
-     Real disp_norm = 0.;
+     //Real disp_norm = 0.;
 
      UInt n = surface.size();
      UInt size = n*n;
 
-   #pragma omp parallel for reduction(+:crit, disp_norm)
+     //#pragma omp parallel for reduction(+:crit, disp_norm)
+#pragma omp parallel for reduction(+:crit)
      for (UInt i = 0; i < size; ++i) {
          crit += std::abs(this->gap(i)*(this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho)));
 
      }
 
 
 
   return crit;
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 
 /* -------------------------------------------------------------------------- */
 
 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;
   const Surface<Real> & gradF = this->functional->getGradF();
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
 
     this->search_direction(i)= gradF(i);
   }
 
   STOPTIMER("updateSearchDirection");
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeOptimalStep() {
   STARTTIMER("computeOptimalStep");
 
-  search_direction.FFTTransform(spectral_search_direction);
+  this->applyInverseInfluenceFunctions(search_direction, projected_direction);
 
   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();
- // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
- // spectral_search_direction -= average;
 
   const 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);
+    denominator += projected_direction(i) * search_direction(i);
   }
 
   Real alpha = numerator / denominator;
 
   STOPTIMER("computeOptimalStep");
   return alpha;
 
 
 
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeTau() {
-  STARTTIMER("computeOptimalStep");
-
-  search_direction.FFTTransform(spectral_search_direction);
-
-  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();
- // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
- // spectral_search_direction -= average;
-
-  const 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;
+  return computeOptimalStep();
 }
 
 /* -------------------------------------------------------------------------- */
 
 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::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 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 BemGigi::computePressures(Real mean_displacement) {
 
   this->computeTractionsFromDisplacements();
 
     UInt n = surface.size();
     UInt size = n*n;
     Real rho = this->functional->getParameter("rho");
     Real surface_energy = this->functional->getParameter("surface_energy");
 
     Real mini=3000;
     for (UInt i = 0; i < size; ++i) {
    if (gap(i)<rho){
     if (this->surface_tractions(i)+surface_energy/rho <mini  )
    mini= this->surface_tractions(i)+surface_energy/rho ;
     }
    else
    { if (this->surface_tractions(i)<mini  )
            mini= this->surface_tractions(i) ;}
   }
 
   this->surface_tractions-=mini;
 
 }
 
 __END_TAMAAS__
diff --git a/src/bem_gigi.hh b/src/bem_gigi.hh
index 03d7edb..da40aa6 100644
--- a/src/bem_gigi.hh
+++ b/src/bem_gigi.hh
@@ -1,98 +1,96 @@
 /**
  *
  * @author Valentine Rey <valentine.rey@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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #ifndef BEM_GIGI_H
 #define BEM_GIGI_H
 /* -------------------------------------------------------------------------- */
 #include "bem_fft_base.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 class BemGigi : public BemFFTBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
 
   BemGigi(Surface<Real> & p) :
     BemFFTBase(p),
     search_direction(p.size(), p.getL()),
-    spectral_search_direction(p.size(), p.getL()),
-     surface_spectral_r(p.size(), p.getL()),
+    projected_direction(p.size(), p.getL()),
     mean_disp_search_direction(p.size(), p.getL()) {
     e_star = 1.;
     max_iterations = 10000;
   };
   virtual ~BemGigi(){};
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   Real computeEquilibrium(Real epsilon, Real mean_displacement);
   void updateSearchDirection(Real factor);
   void updateDisplacements(Real tau);
 
   Real computeOptimalStep();
   void emptyOverlap();
   void enforceMeanDisplacement(Real mean_displacement);
   void computePressures(Real mean_displacement);
   Real computeStoppingCriterion();
   void optimizeToMeanDisplacement(Real imposed_mean);
   Real positiveGapAverage(Real shift);
 
 
   Real computeTau();
 
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
 
   //! get the number of iterations
   UInt getNbIterations() const{return this->nb_iterations;};
   //! get the convergence steps of iterations
   const std::vector<Real> & getConvergenceIterations() const{return this->convergence_iterations;};
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
 
   //! exploration direction for the Gigi algo
   Surface<Real> search_direction;
   //! projected exploration direction for the Gigi algo (fourier version)
-  SurfaceComplex<Real> spectral_search_direction;
-  SurfaceComplex<Real>  surface_spectral_r;
+  Surface<Real> projected_direction;
   //! exploration direction for optimization of mean displacement
   Surface<Real> mean_disp_search_direction;
 };
 
 __END_TAMAAS__
 
 #endif /* BEM_GIGI_H */
diff --git a/src/bem_penalty.cpp b/src/bem_penalty.cpp
index e603bf3..bfbadc1 100644
--- a/src/bem_penalty.cpp
+++ b/src/bem_penalty.cpp
@@ -1,247 +1,233 @@
 /**
  *
  * @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_penalty.hh"
 #include <iostream>
 #include <sstream>
 #include <fstream>
 #include <iomanip>
 #include <sstream>
 #include <cmath>
 /* -------------------------------------------------------------------------- */
 #define TIMER
 #include "surface_timer.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 Real BemPenalty::computeEquilibrium(Real epsilon,
                                  Real mean_displacement, Real penalization_parameter) {
 
   this->computeSpectralInfluenceOverDisplacement();
 
     this->search_direction = 0.;
     this->true_displacements =0;
     this->true_displacements = SurfaceStatistics::computeMaximum(this->surface);
      this->computeGaps();
 
 
 
     Real f = 1e300;
 
   std::cout << "moyenne deplacement "<< SurfaceStatistics::computeAverage(this->true_displacements, 0)<< std::endl;
 
     convergence_iterations.clear();
     nb_iterations = 0;
 
    max_iterations =5000;
 
   while (f> epsilon && nb_iterations++ < max_iterations) {
 
       this->functional->computeGradFU();
     this->computeSearchDirection(mean_displacement,penalization_parameter);
 
     Real alpha=0.0001;
 
        this->old_displacements=this->true_displacements;
     this->updateUnknown( alpha, mean_displacement);
       this->computeGaps();
    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);
   }
 
 
    this->computePressures(mean_displacement);
 
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemPenalty::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->search_direction(i)*this->search_direction(i);
       disp_norm += (true_displacements(i)*true_displacements(i));
     }
 
 
 
  return crit / disp_norm;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemPenalty::computeSearchDirection(Real mean_displacement,Real penalization_parameter) {
   STARTTIMER("computeOptimalStep");
 
 
   UInt n = surface.size();
   UInt size = n*n;
 
  const Surface<Real> & gradF = this->functional->getGradF();
 
 #pragma omp parallel for
   for (UInt i = 1; i < size; ++i) {
   this->search_direction(i)=gradF(i);
 
 
 
       if (gap(i)<0)
       {
          this->search_direction(i)+=penalization_parameter*gap(i);
 
 
 
       }
 }
 
 
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemPenalty::computeOptimalStep() {
   STARTTIMER("computeOptimalStep");
 
-  search_direction.FFTTransform(spectral_search_direction);
-
+  this->applyInverseInfluenceFunctions(search_direction, surface_r);
   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();
- // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
-//spectral_search_direction -= average;
-
-
-
  Real numerator = 0., denominator = 0.;
 
 #pragma omp parallel for reduction(+: numerator, denominator)
  for (UInt i = 0; i < size; ++i) {
    numerator += search_direction(i) * search_direction(i);
-    denominator += spectral_search_direction(i).real() * search_direction(i);
+   denominator += surface_r(i) * search_direction(i);
  }
 
  Real alpha = numerator / denominator;
 
   STOPTIMER("computeOptimalStep");
  return alpha;
 
 
 
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 
 void BemPenalty::updateUnknown(Real alpha, Real mean_displacement) {
   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);
 
   }
 
 Real moyenne=SurfaceStatistics::computeAverage(this->true_displacements, 0);
 
 for (UInt i = 0; i < size; ++i) {
   this->true_displacements(i) =this->true_displacements(i)-moyenne+mean_displacement;
 
 }
 
   STOPTIMER("updateDisplacements");
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemPenalty::computePressures(Real mean_displacement) {
 
   this->computeTractionsFromDisplacements();
 
     UInt n = surface.size();
     UInt size = n*n;
     Real rho = this->functional->getParameter("rho");
     Real surface_energy = this->functional->getParameter("surface_energy");
 
 
    Real mini=3000;
     for (UInt i = 0; i < size; ++i) {
    if (gap(i)<rho){
     if (this->surface_tractions(i)+surface_energy/rho <mini  )
    mini= this->surface_tractions(i)+surface_energy/rho ;
     }
    else
   { if (this->surface_tractions(i)<mini  )
            mini= this->surface_tractions(i) ;}
   //}
 
         //if (this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho)<mini  )
                //  mini= this->surface_tractions(i) +surface_energy/rho * exp(- (gap(i)) / rho);
   }
 
   this->surface_tractions-=mini;
 
 
 
 
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
diff --git a/src/bem_penalty.hh b/src/bem_penalty.hh
index b5f9d0f..b60e6eb 100644
--- a/src/bem_penalty.hh
+++ b/src/bem_penalty.hh
@@ -1,99 +1,93 @@
 /**
  *
  * @author Valentine Rey <valentine.rey@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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #ifndef BEM_PENALTY_H
 #define BEM_PENALTY_H
 /* -------------------------------------------------------------------------- */
 #include "bem_fft_base.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 class BemPenalty : public BemFFTBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
 
   BemPenalty(Surface<Real> & p) :
     BemFFTBase(p),
     search_direction(p.size(), p.getL()),
-    spectral_search_direction(p.size(), p.getL()),
     surface_t(p.size(),p.getL()),
     surface_r(p.size(),p.getL()),
-    surface_spectral_r(p.size(),p.getL()),
     mean_disp_search_direction(p.size(), p.getL()) {
     e_star = 1.;
     max_iterations = 10000;
   };
   virtual ~BemPenalty(){};
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   //! compute the equilibrium situation
   Real computeEquilibrium(Real epsilon,Real pressure,Real penalization_parameter);
   void computeSearchDirection(Real mean_displacement,Real penalization_parameter);
   void updateUnknown(Real tau, Real mean_displacement);
   Real computeStoppingCriterion();
     void computePressures(Real mean_displacement);
     Real computeOptimalStep();
 
 
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
 
   //! get the number of iterations
   UInt getNbIterations() const{return this->nb_iterations;};
   //! get the convergence steps of iterations
   const std::vector<Real> & getConvergenceIterations() const{return this->convergence_iterations;};
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
 
   //! exploration direction for the Gigi algo
   Surface<Real> search_direction;
-  //! projected exploration direction for the Gigi algo (fourier version)
-  SurfaceComplex<Real> spectral_search_direction;
   //! exploration direction for the polonski algo
   Surface<Real> surface_t;
   //! projected exploration direction for the polonski algo
   Surface<Real> surface_r;
-  //! projected exploration direction for the polonski algo (fourier version)
-  SurfaceComplex<Real> surface_spectral_r;
   //! exploration direction for optimization of mean displacement
   Surface<Real> mean_disp_search_direction;
 };
 
 __END_TAMAAS__
 
 #endif /* BEM_PENALTY_H */
diff --git a/src/bem_uzawa.cpp b/src/bem_uzawa.cpp
index 1d9077c..282ff7a 100644
--- a/src/bem_uzawa.cpp
+++ b/src/bem_uzawa.cpp
@@ -1,329 +1,315 @@
 /**
  *
  * @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_uzawa.hh"
 #include <iostream>
 #include <sstream>
 #include <fstream>
 #include <iomanip>
 #include <sstream>
 #include <cmath>
 /* -------------------------------------------------------------------------- */
 #define TIMER
 #include "surface_timer.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 Real BemUzawa::computeEquilibrium(Real epsilon,
 				  Real mean_displacement, Real penalization_parameter) {
 
   this->computeSpectralInfluenceOverDisplacement();
 
   this->search_direction = 0.;
   this->true_displacements = 0;
   this->true_displacements = mean_displacement;
   this->computeGaps();
   this->multipliers = 0;
 
 
   Real f_2 = 1e300;
   Real new_pen = 1e300;
   Real old_pen = 1e300;
 
   convergence_iterations.clear();
   UInt nb_iterations_1 = 0;
   UInt nb_iterations_2 = 0;
 
 
   Real max_iterations_1 = 500;
   Real   max_iterations_2 = 1000;
 
   while (new_pen> epsilon && nb_iterations_1++ < max_iterations_1) {
     std::cout << "iter ext "<< nb_iterations_1 << std::endl;
 
 
     while (f_2> 1e-12 && nb_iterations_2++ < max_iterations_2) {
       std::cout << "iter int "<< nb_iterations_2 << std::endl;
       this->functional->computeGradFU();
       this->computeSearchDirection(mean_displacement,penalization_parameter);
 
       Real  alpha=1/(10*penalization_parameter);
       this->old_displacements=this->true_displacements;
       this->updateUnknown( alpha, mean_displacement);
       this->computeGaps();
 
       f_2 = computeStoppingCriterion();
 
       if (nb_iterations_2 % dump_freq == 0) {
 	std::cout << std::scientific << std::setprecision(10)
 		  << nb_iterations << " "
 		  << f_2 << std::fixed << std::endl;
       }
     }
 
     this->updateMultipliers(penalization_parameter);
 
     old_pen=computeInterpenetration(this->old_displacements);
     new_pen=computeInterpenetration(this->true_displacements);
     std::cout << "old penetration is  "<< old_pen << std::endl;
     std::cout << "new penetration is  "<< new_pen << std::endl;
 
     penalization_parameter=this->updatePenalization(penalization_parameter, old_pen, new_pen);
     //to avoid ill conditioned system
     if (penalization_parameter>1.0e8) new_pen=0.;
     nb_iterations_2 = 0;
     f_2 = 1e300;
     std::cout << "penalization is  "<< penalization_parameter << std::endl;
 
   }
   this->computeGaps();
   this->computePressures();
 
 
 
   return new_pen;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemUzawa::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->search_direction(i))*(this->search_direction(i));
     disp_norm += (true_displacements(i)*true_displacements(i));
   }
 
 
 
   return crit / disp_norm;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemUzawa::computeSearchDirection(Real mean_displacement,Real penalization_parameter) {
   STARTTIMER("computeOptimalStep");
 
 
   UInt n = surface.size();
   UInt size = n*n;
 
   const Surface<Real> & gradF = this->functional->getGradF();
 
 #pragma omp parallel for
   for (UInt i = 1; i < size; ++i) {
     this->search_direction(i)=gradF(i);
 
     if (this->multipliers(i)+penalization_parameter*gap(i)<=0){
 
 
       this->search_direction(i)=this->search_direction(i)+  this->multipliers(i)+penalization_parameter*gap(i);
     }
 
   }
 
 
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemUzawa::computeOptimalStep() {
   STARTTIMER("computeOptimalStep");
 
-  search_direction.FFTTransform(spectral_search_direction);
-
+  this->applyInverseInfluenceFunctions(search_direction, surface_r);
   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();
-  // Real average = SurfaceStatistics::computeAverage(spectral_search_direction.real(), 0);
-  //spectral_search_direction -= average;
-
-
-
   Real numerator = 0., denominator = 0.;
 
 #pragma omp parallel for reduction(+: numerator, denominator)
   for (UInt i = 0; i < size; ++i) {
     numerator += search_direction(i) * search_direction(i);
-    denominator += spectral_search_direction(i).real() * search_direction(i);
+    denominator += surface_r(i) * search_direction(i);
   }
 
   Real alpha = numerator / denominator;
 
   STOPTIMER("computeOptimalStep");
   return alpha;
 
 
 
 }
 
 
 /* -------------------------------------------------------------------------- */
 void BemUzawa::updateMultipliers(Real penalization_parameter) {
   STARTTIMER("updateMultipliers");
 
   UInt n = surface.size();
   UInt size = n*n;
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     if (this->multipliers(i)+penalization_parameter*gap(i)>0.){
       this->multipliers(i) = 0;
     }
     else{
       this->multipliers(i)+=penalization_parameter*gap(i);
     }
 
 
   }
 
 
 
   STOPTIMER("updateMultipliers");
 }
 
 /* -------------------------------------------------------------------------- */
 Real BemUzawa::updatePenalization(Real penalization_parameter, Real old_pen, Real new_pen) {
   STARTTIMER("updatePenalization");
 
 
   if (new_pen>0.25*old_pen){
     penalization_parameter*=5;
   }
 
 
   return penalization_parameter;
 
 
 
 
   STOPTIMER("updatePenalization");
 }
 
 
 /* -------------------------------------------------------------------------- */
 Real BemUzawa::computeInterpenetration( Surface<Real> & displacements ) {
   STARTTIMER("updateMultipliers");
 
   UInt n = surface.size();
   UInt size = n*n;
   Real res=0.;
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     if (0>this->true_displacements(i)-surface(i)){
       res+=(this->true_displacements(i)-surface(i))*(this->true_displacements(i)-surface(i));
     }
 
   }
 
   return sqrt(res);
 
 
 
   STOPTIMER("updateMultipliers");
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 
 void BemUzawa::updateUnknown(Real alpha, Real mean_displacement) {
   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);
 
   }
 
   Real moyenne=SurfaceStatistics::computeAverage(this->true_displacements, 0);
 
   for (UInt i = 0; i < size; ++i) {
     this->true_displacements(i) =this->true_displacements(i)-moyenne+mean_displacement;
 
   }
 
 
 
   STOPTIMER("updateDisplacements");
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 void BemUzawa::computePressures() {
   this->applyInverseInfluenceFunctions(this->true_displacements, this->surface_tractions);
 
 
     UInt n = surface.size();
     UInt size = n*n;
     Real rho = this->functional->getParameter("rho");
     Real surface_energy = this->functional->getParameter("surface_energy");
 
 
    Real mini=3000;
     for (UInt i = 0; i < size; ++i) {
    if (gap(i)<rho){
     if (this->surface_tractions(i)+surface_energy/rho <mini  )
    mini= this->surface_tractions(i)+surface_energy/rho ;
     }
    else
   { if (this->surface_tractions(i)<mini  )
            mini= this->surface_tractions(i) ;}
   //}
 
         //if (this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho)<mini  )
                //  mini= this->surface_tractions(i) +surface_energy/rho * exp(- (gap(i)) / rho);
   }
 
   this->surface_tractions-=mini;
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
diff --git a/src/bem_uzawa.hh b/src/bem_uzawa.hh
index 2ec22d2..2478c32 100644
--- a/src/bem_uzawa.hh
+++ b/src/bem_uzawa.hh
@@ -1,100 +1,94 @@
 /**
  *
  * @author Valentine Rey <valentine.rey@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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #ifndef BEM_UZAWA_H
 #define BEM_UZAWA_H
 /* -------------------------------------------------------------------------- */
 #include "bem_fft_base.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 class BemUzawa : public BemFFTBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
 
   BemUzawa(Surface<Real> & p) :
     BemFFTBase(p),
     search_direction(p.size(), p.getL()),
-    spectral_search_direction(p.size(), p.getL()),
     multipliers(p.size(),p.getL()),
     surface_r(p.size(),p.getL()),
-    surface_spectral_r(p.size(),p.getL()),
     mean_disp_search_direction(p.size(), p.getL()) {
     e_star = 1.;
     max_iterations = 10000;
   };
   virtual ~BemUzawa(){};
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   //! compute the equilibrium situation
   Real computeEquilibrium(Real epsilon,Real pressure,Real penalization_parameter);
   void computeSearchDirection(Real mean_displacement,Real penalization_parameter);
   void updateMultipliers(Real penalization_parameter);
   Real updatePenalization(Real penalization_parameter, Real old_pen, Real new_pen);
   Real computeInterpenetration( Surface<Real> & displacements ) ;
   void updateUnknown(Real tau, Real mean_displacement);
   Real computeStoppingCriterion();
   void computePressures();
   Real computeOptimalStep();
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
 
   //! get the number of iterations
   UInt getNbIterations() const{return this->nb_iterations;};
   //! get the convergence steps of iterations
   const std::vector<Real> & getConvergenceIterations() const{return this->convergence_iterations;};
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
 
   //! exploration direction for the Gigi algo
   Surface<Real> search_direction;
-  //! projected exploration direction for the Gigi algo (fourier version)
-  SurfaceComplex<Real> spectral_search_direction;
   //! exploration direction for the polonski algo
   Surface<Real> multipliers;
   //! projected exploration direction for the polonski algo
   Surface<Real> surface_r;
-  //! projected exploration direction for the polonski algo (fourier version)
-  SurfaceComplex<Real> surface_spectral_r;
   //! exploration direction for optimization of mean displacement
   Surface<Real> mean_disp_search_direction;
 };
 
 __END_TAMAAS__
 
 #endif /* BEM_Uzawa_H */