diff --git a/src/bem_gigi.cpp b/src/bem_gigi.cpp
index c0ae58f..4c683b5 100644
--- a/src/bem_gigi.cpp
+++ b/src/bem_gigi.cpp
@@ -1,478 +1,480 @@
 /**
  *
  * @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);
 
 
 
   std::cout << "moyenne deplacement "<< SurfaceStatistics::computeAverage(this->true_displacements, 0)<< std::endl;
 
 
   convergence_iterations.clear();
   nb_iterations = 0;
  UInt crit2=10;
 
 
    max_iterations =2000;
 
 
   while (f> epsilon && nb_iterations++ < max_iterations) {
 
     this->computeGaps();
 
 
     this->functional->computeGradF();
 
     Real R = this->functional->computeGradientNorm();
 
     std::cout << "norme gradient "<< R << std::endl;
 
 
 
     this->updateSearchDirection(R / Rold);
 
     Rold = R;
     Real alpha = this->computeOptimalStep();
     //alpha=0.001;
         std::cout << "alpha vaut "<< alpha<< std::endl;
 
     this->old_displacements=this->true_displacements;
     this->updateDisplacements(alpha);
 
   //std::cout << "moyenne deplacement "<< SurfaceStatistics::computeAverage(this->true_displacements, 0)<< std::endl;
 
     this->computeGaps();
 
 
 
     this->optimizeToMeanDisplacement(mean_displacement); //espace admissible
     this->computeGaps();
 
 
 
   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 ) {
       ++inner_loop_cpt;
        this->computeGaps();
       this->functional->computeGradF();
       this->updateSearchDirection(R / Rold);
 
       Real tau = this->computeTau();
-      Gold = G;
+      //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);
 
   }
 
 
 
    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)= 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))-0.002/2.071e-3 * exp(- gap(i) / 2.071e-3)*gap(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;
 
 
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 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;
 
 
       }
   }
 
 
 
 
 }
diff --git a/src/bem_penalty.hh b/src/bem_penalty.hh
index b82e396..e1bf7a9 100644
--- a/src/bem_penalty.hh
+++ b/src/bem_penalty.hh
@@ -1,95 +1,95 @@
 /**
  *
  * @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"
 /* -------------------------------------------------------------------------- */
 
 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()),
-    search_direction(p.size(), p.getL()),
-    spectral_search_direction(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;
 };
 
 #endif /* BEM_PENALTY_H */
diff --git a/src/surface.cpp b/src/surface.cpp
index d715713..2a99b2f 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -1,449 +1,448 @@
 /**
  *
  * @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 <iostream>
 #include <fstream>
 #include <sstream>
 #include <string>
 #include <cmath>
 #include <random>
 /* -------------------------------------------------------------------------- */
 #include "surface.hh"
 #include "surface_statistics.hh"
 /* -------------------------------------------------------------------------- */
 #ifdef USING_IOHELPER
 #include "dumper_paraview.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 std::string fname = "surface_height.pdf";
 /* -------------------------------------------------------------------------- */
 
 
 
 template <typename T>
 void Surface<T>::expandByLinearInterpolation() {
   UInt n = this->size();
   Surface * old_s = new Surface(n,1.);
   int new_n = n*2;
   *old_s = *this; 
   this->setGridSize(new_n);
   for (UInt i=0; i<n; i++) {
     for (UInt j=0; j<n; j++) {
       (*this)(i,j) = (*old_s)(i,j);
     }
   }
   delete old_s;
 }
 
 Real sinc(Real x){
   return sin(M_PI*x)/(M_PI*x);
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByShannonInterpolation(int factor) {
 
   this->FFTTransform();
   UInt n = this->size();
   UInt N = std::pow(2,factor)*this->size();
   std::cout << N << " " << factor << std::endl;
   Surface s = Surface(N,1.);
 
 
   for (UInt i = 0; i < n/2; ++i) 
     for (UInt j = 0; j < n/2; ++j) {
       s(i,j) = this->at(i,j);
     }
 
   for (UInt i = n/2; i < n; ++i) 
     for (UInt j = 0; j < n/2; ++j) {
       s(N-n+i,j) = this->at(i,j);
     }
 
   for (UInt i = n/2; i < n; ++i) 
     for (UInt j = n/2; j < n; ++j) {
       s(N-n+i,N-n+j) = this->at(i,j);
     }
 
   for (UInt i = 0; i < n/2; ++i) 
     for (UInt j = n/2; j < n; ++j) {
       s(i,N-n+j) = this->at(i,j);
     }
 
   //  s.dumpToTextFile<false,true>("temp.csv");
 
   s.FFTITransform();
 
   this->setGridSize(N);
   Real f = std::pow(2,factor);
   f = f*f;
 
   for (UInt i = 0; i < N*N; ++i) 
     this->at(i) = f*s(i);
   
  
 //   int new_n = 2*n;
 //   Surface s(*this);
 //   this->setGridSize(new_n);
 
 //   // put the original signal with zeros for to be interpolated points
 //   for (int i=0; i<n; i++) 
 //     for (int j=0; j<n; j++) 
 //       this->at(2*i,j*2) = s(i,j);
 
 //   // do shannon interpolation
 //   int rcut = 64;
 
 //   //interpolate lines
 //   for (int i=0; i<new_n; i++) {
 //     for (int j=0; j<new_n; j++) {
 //       if (i%2 == 0) continue;
 //       // Real loop on local neighbors to interpolate
 //       for (int k=-rcut; k<rcut+1; k++) {
 // 	int i2 = i + k;
 // 	if (i2 < 0) i2 += new_n;
 // 	if (i2 >= new_n) i2 -= new_n;
 // 	if (i2%2) continue;
 	
 // 	Real d = fabs(k);
 // 	s[i+j*new_n] += s(i2,j)*sinc(d/2);
 //       }
 //     }
 //   }
 
 //   //interpolate columns
 //   for (int i=0; i<new_n; i++) {
 //     for (int j=0; j<new_n; j++) {
 //       if (j%2 == 0) continue;
 //       // Real loop on local neighbors to interpolate
 //       for (int k=-rcut; k<rcut+1; k++) {
 // 	int j2 = j + k;
 // 	if (j2 < 0) j2 += new_n;
 // 	if (j2 >= new_n) j2 -= new_n;
 // 	if (j2%2) continue;
 	
 // 	Real d = fabs(k);
 // 	s[i+j*new_n] += complex(s[i+j2*new_n].real()*sinc(d/2),
 // 				s[i+j2*new_n].imag()*sinc(d/2));
 //       }
 //     }
 //   }
 
   
 // //   for (int i=0; i<n; i++) {
 // //     for (int j=0; j<n; j++) {
 // //       s[2*i+j*2*new_n].real() = 0;
 // //       s[2*i+j*2*new_n].im = 0;
 // //     }
 // //   }
 
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByMirroring() {
   UInt n = this->size();
   UInt new_n = n*2;
   Surface s(*this);
   this->setGridSize(new_n);
 
   for (UInt i=0; i<new_n; i++) {
     for (UInt j=0; j<new_n; j++) {
       UInt i1 = i;
       if (i1 > n) i1 = new_n - i1;
       UInt j1 = j;
       if (j1 > n) j1 = new_n - j1;
       if (i == n || j == n)  {
 	this->at(i,j) = 0.;
 	continue;
       }
       this->at(i,j) = s(i1,j1);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByPuttingZeros() {
   UInt n = this->size();
   UInt new_n = n*2;
   Surface s(*this);
   this->setGridSize(new_n);
 
   for (UInt i=0; i< n; i++) {
     for (UInt j=0; j< n; j++) {
       this->at(2*i,j*2) = s(i,j);
     }
   }
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::contractByTakingSubRegion() {
   UInt n = this->size();
   Surface s(*this);
   this->setGridSize(n/2);
 
   for (UInt i=0; i < n/2; i++) {
     for (UInt j=0; j < n/2; j++) {
       this->at(i,j) = s(i,j);
     }
   }
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::contractByIgnoringMidPoints() {
   UInt n = this->size();
   Surface s(*this);
   this->setGridSize(n/2);
 
   for (UInt i=0; i<n; i++) {
     for (UInt j=0; j<n; j++) {
       if (i%2) continue;
       if (j%2) continue;
       this->at(i/2,j/2) = s(i,j);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 // void Surface<T>::dumpToTextFile(string filename, bool pr) {
 
 //   cout << "Writing surface file " << filename << endl;
 //   ofstream file(filename.c_str());
 //   file.precision(15);
 
 //   if (pr)
 //   {
 //     for (int i=0; i<n; i++) {
 //       for (int j=0; j<n; j++) {
 //         file << i << " " << j << " ";
 //         if (heights[i+j*n].real() < 1e-10)
 //           file << scientific << -1e10 << " " << 0;
 //         else
 //           file << scientific << heights[i+j*n].real() << " " << heights[i+j*n].im;
 //         file << endl;
 //       }
 //       file << endl;
 //     }
 //   }
 //   else
 //   {
 //     for (int i=0; i<n; i++) {
 //       for (int j=0; j<n; j++) {
 //         file << i << " " << j << " ";
 //         file << scientific << heights[i+j*n].real() << " " << heights[i+j*n].im;
 //         file << endl;
 //       }
 //     }
 //   }
 //   file.close();
   
 // }
 
 /* -------------------------------------------------------------------------- */
 
 
 
 // Surface * Surface<T>::loadFromTextFile(std::string filename) {
 
   // std::string line;
   // std::ifstream myfile (filename.c_str());
 
   // if (!myfile.is_open()){
   //   SURFACE_FATAL("cannot open file " << filename);
   // }
 
   // Surface * my_surface = NULL;
   // UInt n = 0;
   // if (myfile.is_open())
   //   {
   //     while (! myfile.eof() )
   // 	{
   // 	  getline (myfile,line);
   // 	  if (line[0] == '#')
   // 	    continue;
   // 	  if (line.size() == 0)
   // 	    continue;
 
   // 	  std::stringstream s;
   // 	  s << line;
   // 	  int i,j;
   // 	  s >> i;
   // 	  s >> j;
   // 	  if (n < i) n = i;
   // 	  if (n < j) n = j;
   // 	}
   //     ++n;
   //     std::cout << "read " << n << " x " << n << " grid points" << std::endl;
       
   //     myfile.clear();              // forget we hit the end of file
   //     myfile.seekg(0, std::ios::beg);   // move to the start of the file
 
   //     my_surface = new Surface(n);
 
   //     while (! myfile.eof() )
   // 	{
   // 	  getline (myfile,line);
   // 	  if (line[0] == '#')
   // 	    continue;
   // 	  if (line.size() == 0)
   // 	    continue;
 
   // 	  std::stringstream s;
   // 	  s << line;
   // 	  int i,j;
   // 	  Real h,hIm;
   // 	  s >> i;
   // 	  s >> j;
   // 	  s >> h;
   // 	  s >> hIm;
   // 	  my_surface->heights[i+j*n] = complex(h,hIm);
   // 	}
   //     myfile.close();
   //   }
 //   return my_surface;
 // }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::dumpToParaview(std::string filename) const {
 #ifdef USING_IOHELPER
 
   UInt n = this->size();
   const Surface & data = *this; 
   /* build 3D coordinates array to be usable with iohelper */
   //  cerr << "allocate coords for paraview" << endl;
   Real *  coords = new Real[n*n*3];
   //Real *  coords = new Real[(2*n-1)*(2*n-1)*3];
   //cerr << "allocate zcoords for paraview" << endl;
   Real *  zcoords = new Real[n*n*3];
   //Real *  zcoords = new Real[(2*n-1)*(2*n-1)*3];
   //cerr << "allocate zcoords2 for paraview" << endl;
   Real *  zcoords2 = new Real[n*n*3];
   //  Real *  zcoords2 = new Real[(2*n-1)*(2*n-1)*3];
   int index = 0;
   // for (int i=-n+1; i<n; i++) {
   //   for (int j=-n+1; j<n; j++) {
   for (int i=0; i<n; i++) {
     for (int j=0; j<n; j++) {
       coords[3*index] = 1.*i/n;
       coords[3*index+1] = 1.*j/n;
       coords[3*index+2] = 0;
       zcoords[3*index] = 0;
       zcoords[3*index+1] = 0;
       zcoords2[3*index] = 0;
       zcoords2[3*index+1] = 0;
       int i1 = (i%n);
       int j1 = (j%n);
       if (i1 < 0) i1 += n;
       if (j1 < 0) j1 += n;
       zcoords[3*index+2] = data(i1,j1).real();
       zcoords2[3*index+2] = data(i1,j1).imag();
       ++index;
     }
   }
   
   iohelper::DumperParaview dumper;
   dumper.setPoints(coords, 3, index, filename.c_str());
   dumper.setConnectivity(NULL, iohelper::POINT_SET, 1, iohelper::FORTRAN_MODE);
   dumper.addNodeDataField("data.re",zcoords,3,index);
   dumper.addNodeDataField("data.im",zcoords2,3,index);
   dumper.setPrefix("./paraview");
   dumper.init();
   dumper.dump();
 
   delete [] coords;
   delete [] zcoords;
   delete [] zcoords2;
 #endif //USING_IOHELPER
 }
 /* -------------------------------------------------------------------------- */
 
 
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 Surface<T>::Surface(UInt a, Real L) : Map2dSquare<T>(a,L){
   this->setGridSize(a);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 Surface<T>::~Surface(){
 }
 
 /* -------------------------------------------------------------------------- */
 
 // void Surface<T>::copy(Surface & s){
 //   this->n = s.n;
 //   this->L = s.L;
 //   this->heights = s.heights;
 // }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::recenterHeights(){
   UInt n = this->size();
   Real zmean = SurfaceStatistics::computeAverage(*this);
   
   for (UInt i = 0 ; i < n*n ; ++i)
     this->at(i) = this->at(i) - zmean;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::generateWhiteNoiseSurface(Surface<Real> & sources, long random_seed){
 
     //set the random gaussian generator
       std::mt19937 gen(random_seed);
       std::normal_distribution<> gaussian_generator;
       Surface<Real> * origins;
        origins = new Surface<Real>(1024,1.);
       UInt n = sources.size();
 
       for (UInt i = 0 ; i < 1024 ; ++i)
         for (UInt j = 0 ; j < 1024 ; ++j){
           origins->at(i,j) = gaussian_generator(gen);
         }
 
       for (UInt i = 0 ; i < n ; ++i)
         for (UInt j = 0 ; j < n ; ++j){
           sources(i,j) = origins->at(i*1024/n,j*1024/n);
         }
 
       delete origins;
-      std::cout << "\t>> coucou from white noise" << std::endl;
 }
 /* -------------------------------------------------------------------------- */
 
 template class Surface<Real>;
 template class Surface<unsigned long>;
 template class Surface<UInt>;
 template class Surface<int>;
 template class Surface<bool>;