diff --git a/python/bem.i b/python/bem.i
index e8d0b12..eb21902 100644
--- a/python/bem.i
+++ b/python/bem.i
@@ -1,57 +1,58 @@
 /**
  *
  * @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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 
 
 /* -------------------------------------------------------------------------- */
 /* Wrapping bem capabilities
 /* -------------------------------------------------------------------------- */
 
 %include "bem_interface.hh"
 %include "bem_fft_base.hh"
 %include "bem_kato.hh"
 %include "bem_polonski.hh"
 %include "bem_gigi.hh"
 %include "bem_penalty.hh"
 %include "bem_uzawa.hh"
 %include "bem_gigipol.hh"
 %include "bem_grid.hh"
 %include "bem_grid_kato.hh"
 %include "bem_grid_polonski.hh"
 %include "bem_grid_teboulle.hh"
 %include "functional.hh"
 %include "meta_functional.hh"
 %include "elastic_energy_functional.hh"
 %include "complimentary_term_functional.hh"
 %include "exponential_adhesion_functional.hh"
+%include "squared_exponential_adhesion_functional.hh"
 %include "maugis_adhesion_functional.hh"
 %include "fftransform.hh"
 %include "fftransform_fftw.hh"
 
 namespace tamaas {
   %template(FFTransformReal) FFTransform<Real,2>;
   %template(FFTransformFFTWReal) FFTransformFFTW<Real,2>;
 }
diff --git a/python/surface.i b/python/surface.i
index 68a24ca..683cefb 100644
--- a/python/surface.i
+++ b/python/surface.i
@@ -1,379 +1,380 @@
 /**
  *
  * @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 <complex>
 #include "map_2d.hh"
 #include "surface.hh"
 #include "types.hh"
 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 #include <numpy/arrayobject.h>
 #include "surface_generator_filter_fft.hh"
 #include "smart_pointer.hh"
 #include "surface_statistics.hh"
 #include "bem_kato.hh"
 #include "bem_polonski.hh"
 #include "bem_penalty.hh"
 #include "bem_gigi.hh"
 #include "bem_gigipol.hh"
 #include "bem_uzawa.hh"
 #include "bem_fft_base.hh"
 #include "bem_grid.hh"
 #include "bem_grid_polonski.hh"
 #include "bem_grid_kato.hh"
 #include "bem_grid_teboulle.hh"
 #include "contact_area.hh"
 #include "contact_cluster_collection.hh"
 #include "functional.hh"
 #include "meta_functional.hh"
 #include "elastic_energy_functional.hh"
 #include "complimentary_term_functional.hh"
 #include "exponential_adhesion_functional.hh"
+#include "squared_exponential_adhesion_functional.hh"
 #include "maugis_adhesion_functional.hh"
 #include "surface_timer.hh"
 #include "fftransform_fftw.hh"
 %}
 %include "cpointer.i"
 %include "typemaps.i"
 %include "numpy.i"
 %init %{
   import_array();
 %}
 
  /* -------------------------------------------------------------------------- */
 /* Wrapping returned pointers
 /* -------------------------------------------------------------------------- */
 //%include "typemaps.i"
 %include "std_string.i"
 %include "std_vector.i"
 %include "smart_pointer.hh"
 %template(VectorReal) std::vector<double>;
 //%pointer_functions(int, intp);
 //%pointer_functions(long, longp);
 //%pointer_functions(unsigned int, uintp);
 //%pointer_functions(double, doublep);
 //%pointer_functions(float, floatp);
 //%pointer_functions(std::complex<double>, comlpexp);
 
 %define %my_pointer_class(NAME,TYPE)
 
 %template(SmartPointer##NAME) SmartPointer< TYPE >;
 
 %typemap(out) TYPE & {
 #define TOREPLACE SWIGTYPE_p_SmartPointerT_##NAME##_t
   SmartPointer<TYPE> * ptr = new SmartPointer<TYPE>($1);
   $result = SWIG_NewPointerObj(SWIG_as_voidptr(ptr), TOREPLACE, 0 |  0 );
 #undef TOREPLACE
  }
 %enddef
 
 %my_pointer_class(int,int)
 %my_pointer_class(double,double)
 %my_pointer_class(long,long)
 %my_pointer_class(vectorReal,std::vector<tamaas::Real>)
 
 %apply (unsigned int IN_ARRAY1[ANY]) {(tamaas::UInt n[ANY])};
 
 /* -------------------------------------------------------------------------- */
 /* Wrapping Surface class                                                     */
 /* -------------------------------------------------------------------------- */
 
 %include "tamaas.hh"
 
 %ignore tamaas::Grid::iterator::operator++;
 %ignore *::operator=;
 
 %include "grid.hh"
 %include "grid_hermitian.hh"
 namespace tamaas {
   %template (Grid1dReal) Grid<Real, 1>;
   %template (Grid2dReal) Grid<Real, 2>;
   %template (Grid2dComplex) Grid<complex, 2>;
   %template (Grid2dInt) Grid<Int, 2>;
   %template (Grid2dUInt) Grid<UInt, 2>;
   %template (GridHermitian2dReal) GridHermitian<Real, 2>;
 }
 
 %include "types.hh"
 
 %include "grid_hermitian.hh"
 
 %include "map_2d.hh"
 namespace tamaas {
 %template(Map2dComplex) Map2d<complex>;
 %template(Map2dInt)     Map2d<int>;
 %template(Map2dReal)    Map2d<Real>;
 }
 
 %include "map_2d_square.hh"
 namespace tamaas {
 %template(Map2dSquareComplex) Map2dSquare<complex>;
 %template(Map2dSquareInt)     Map2dSquare<int>;
 %template(Map2dSquareReal)    Map2dSquare<Real>;
 }
 
 %include "surface.hh"
 %include "surface_complex.hh"
 
 namespace tamaas {
 %template(SurfaceReal) Surface<Real>;
 %template(SurfaceInt)  Surface<int>;
 %template(SurfaceRealComplex) SurfaceComplex<Real>;
 }
 
 %inline %{
   __BEGIN_TAMAAS__
     template <int required_type>
     void sanityCheck(PyObject * input) {
       if (PyArray_TYPE((PyArrayObject*)input) != required_type)
 	throw Exception("Sanity check: python object stores wrong type");
     }
   __END_TAMAAS__
 %}
 
 %inline %{
   __BEGIN_TAMAAS__
 
   template <typename T, UInt d>
   class GridForPython : public Grid<T, d> {
     public:
     GridForPython(T * wrapped_memory, UInt n[d], Real L[d], UInt nb_components):
     Grid<T, d>() {
       this->__init__(wrapped_memory, n, L);
     }
 
     void __init__(T * w, UInt n[d], UInt nb_components) {
       memcpy(this->n, n, d * sizeof(UInt));
       memset(this->L, 1, d * sizeof(Real));
       this->nb_components = nb_components;
       this->data.wrapMemory(w, this->computeSize());
     }
 
     GridForPython(PyObject * input): Grid<T, d>() {
       if (!PyArray_Check(input)) {
 	PyObject* obj_type = PyObject_Type(input);
 	std::string type = PyString_AsString(obj_type);
 	throw Exception("GridForPython: incompatible input which is not a numpy: " + type);
       }
 
       UInt _n = PyArray_NDIM((PyArrayObject*)input); // numpy number of dimensions
       npy_intp * ndims = PyArray_DIMS((PyArrayObject*)input); // numpy dimensions
       UInt grid_dims[d] = {0}; // tamaas dims
       if (_n >= d) {
 	for (UInt i = 0 ; i < d ; ++i)
 	  grid_dims[i] = ndims[i];
       } else {
 	throw Exception("GridForPython: incompatible input size for numpy array: " + _n);
       }
 
       if (_n > d) this->nb_components = ndims[d];
       else this->nb_components = 1;
 
       PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew(input);
       if (iter == NULL) throw;
       this->__init__(( T *)(iter->dataptr), grid_dims, this->nb_components);
       Py_DECREF(iter);
     }
 
     virtual ~GridForPython() {}
   };
 
   __END_TAMAAS__
 %}
 
 %inline %{
   __BEGIN_TAMAAS__
 template <typename T>
   class SurfaceForPython : public Surface<T>{
 
  public:
 
  SurfaceForPython(T * wrapped_memory,
 		   UInt n,
 		   Real L)
    : Surface<T>(0,L){
     this->__init__(wrapped_memory,n,L);
   };
 
 
   void __init__(T * wrapped_memory,
 		UInt n,
 		Real L){
     this->n[0] = n;
     this->n[1] = n;
     this->data.wrapMemory(wrapped_memory,n*n);
   }
 
   SurfaceForPython(PyObject * input, Real L = 1.)
     : Surface<T>(0,L){
 
     if (!PyArray_Check(input)){
       PyObject* obj_type = PyObject_Type(input);
       std::string type = PyString_AsString(obj_type);
       throw Exception("SurfaceForPython: incompatible input which is not a numpy: " + type);
     }
 
     UInt _n = PyArray_NDIM((PyArrayObject*)input);
     if (_n != 2) SURFACE_FATAL("SurfaceForPython: incompatible numpy dimension " << _n);
     npy_intp * ndims = PyArray_DIMS((PyArrayObject*)input);
     UInt sz = ndims[0];
     UInt sz2 = ndims[1];
     if (sz != sz2) SURFACE_FATAL("SurfaceForPython: incompatible numpy dimensions " << sz << " " << sz2);
     PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew(input);
     if (iter == NULL) throw;
     this->__init__(( T *)(iter->dataptr),sz,L);
     Py_DECREF(iter);
   };
 
   ~SurfaceForPython(){
   };
 
   void resize(UInt new_size){
     if (this->size() != new_size)
       throw Exception("SurfaceForPython: cannot resize a temporary vector");
   }
  };
   __END_TAMAAS__
 %}
 
 %define %grid_typemaps(GRID_CLASS, DATA_TYPE, DATA_TYPECODE, DIM)
 %typemap(out, fragment="NumPy_Fragments") (GRID_CLASS &)
 {
   using namespace tamaas;
   UInt additional_dim = ($1->getNbComponents() == 1)? 0:1;
   UInt dim = DIM+additional_dim;
   npy_intp dims[dim];
   for (UInt i = 0 ; i < DIM ; i++) dims[i] = (npy_intp)$1->sizes()[i];
   if (additional_dim)
     dims[dim-1] = $1->getNbComponents();
 
   DATA_TYPE * data = const_cast<DATA_TYPE *>($1->getInternalData());
   PyObject* obj = PyArray_SimpleNewFromData(dim,
 					    dims,
 					    DATA_TYPECODE,
 					    &data[0]);
   PyArrayObject* array = (PyArrayObject*) obj;
   if (!array) SWIG_fail;
   $result = SWIG_Python_AppendOutput($result, obj);
 }
 
 %typemap(in) tamaas::Grid<DATA_TYPE, DIM> & {
   using namespace tamaas;
   if (!PyArray_Check($input)) {
     PyObject* objectsRepresentation = PyObject_Repr(PyObject_Type($input));
     const char* s = PyString_AsString(objectsRepresentation);
     std::string __type = s;
 
     if ("<class 'tamaas.GridForPython" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else if ("<class 'tamaas.Grid" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else throw tamaas::Exception("typemap: incompatible input which is not a proper type: " + __type + " (DATATYPE = " + std::string("DATA_TYPE") + ")");
   } else {
     $1 = new GridForPython< DATA_TYPE, DIM >($input);
     sanityCheck< DATA_TYPECODE >($input);
   }
 }
 
 %inline %{
   namespace tamaas {
   GRID_CLASS & convertGrid(GRID_CLASS & s){
     return s;
   }
   }
 %}
 %enddef
 
 
 %define %surface_typemaps(SURFACE_CLASS, DATA_TYPE, DATA_TYPECODE)
 
 %typemap(out, fragment="NumPy_Fragments") (SURFACE_CLASS &)
 {
   using namespace tamaas;
   UInt nz = $1->size();
   npy_intp dims[2] = {nz,nz};
   DATA_TYPE * data = const_cast<DATA_TYPE *>($1->getInternalData());
   PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, &data[0]);
   PyArrayObject* array = (PyArrayObject*) obj;
   if (!array) SWIG_fail;
   $result = SWIG_Python_AppendOutput($result, obj);
 }
 
 %typemap(in) tamaas::Surface< DATA_TYPE > & {
   using namespace tamaas;
   if (!PyArray_Check($input)){
     PyObject* objectsRepresentation = PyObject_Repr(PyObject_Type($input));
     const char* s = PyString_AsString(objectsRepresentation);
     std::string __type = s;
 
     if ("<class 'tamaas.SurfaceForPython" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else if ("<class 'tamaas.Surface" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else if ("<class 'tamaas.SurfaceComplex" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else throw tamaas::Exception("typemap: incompatible input which is not a proper type: " + __type + " (DATATYPE = " + std::string("DATA_TYPE") + ")");
   }
   else{
     $1 = new SurfaceForPython< DATA_TYPE >($input);
     sanityCheck< DATA_TYPECODE >($input);
   }
  }
 
 %inline %{
   namespace tamaas {
   SURFACE_CLASS & convertSurface(SURFACE_CLASS & s){
     return s;
   }
   }
 %}
 
 
 %enddef
 
 %surface_typemaps(tamaas::Surface<Real>, Real,       NPY_DOUBLE)
 %surface_typemaps(tamaas::SurfaceComplex<Real>, std::complex<Real>,        NPY_COMPLEX128 )
 %surface_typemaps(tamaas::Surface<unsigned long>, unsigned long, NPY_UINT  )
 %surface_typemaps(tamaas::Surface<UInt>, UInt, NPY_UINT  )
 %surface_typemaps(tamaas::Surface<int>, int, NPY_INT  )
 
 %grid_typemaps(tamaas::Grid1dReal, Real, NPY_DOUBLE, 1)
 %grid_typemaps(tamaas::Grid2dReal, Real, NPY_DOUBLE, 2)
 %grid_typemaps(tamaas::GridHermitian2dReal, tamaas::complex, NPY_COMPLEX128, 2)
 %grid_typemaps(tamaas::Grid2dUInt, UInt, NPY_UINT, 2)
 %grid_typemaps(tamaas::SurfaceComplex<Real>, tamaas::complex, NPY_COMPLEX128, 2)
 
 namespace tamaas {
 %template(SurfaceForPythonReal) SurfaceForPython< Real >;
 }
 
 %include "surface_timer.hh"
diff --git a/src/SConscript b/src/SConscript
index ed1507a..6d197cd 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,77 +1,78 @@
 Import('main_env')
 
 import os
 
 env = main_env
 
 # Lib roughcontact
 generator_list = """
 surface_generator.cpp
 surface_generator_voss.cpp
 surface_generator_crenel.cpp
 surface_generator_ellipsoid.cpp
 surface_generator_filter.cpp
 surface_generator_filter_bessel.cpp
 surface_generator_filter_fft.cpp
 """.split()
 #env.SharedLibrary('Generator',generator_list)
 
 # Lib SURFACE
 surface_list = """
 map_2d.cpp
 map_2d_square.cpp
 surface.cpp
 surface_timer.cpp
 tamaas.cpp
 """.split()
 #env.SharedLibrary('Surface',surface_list)
 
 # Lib PERCOLATION
 percolation_list = """
 cluster_grow.cpp
 contact_area.cpp
 contact_cluster.cpp
 contact_cluster_collection.cpp
 """.split()
 #env.SharedLibrary('Percolation',percolation_list)
 
 # BEM PERCOLATION
 bem_list = """
 bem_kato.cpp
 bem_polonski.cpp
 bem_gigi.cpp
 bem_gigipol.cpp
 bem_penalty.cpp
 bem_uzawa.cpp
 bem_fft_base.cpp
 functional.cpp
 meta_functional.cpp
 elastic_energy_functional.cpp
 exponential_adhesion_functional.cpp
+squared_exponential_adhesion_functional.cpp
 maugis_adhesion_functional.cpp
 complimentary_term_functional.cpp
 fftransform.cpp
 fftransform_fftw.cpp
 fft_plan_manager.cpp
 grid.cpp
 grid_hermitian.cpp
 types.cpp
 bem_grid.cpp
 bem_grid_polonski.cpp
 bem_grid_kato.cpp
 bem_grid_teboulle.cpp
 """.split()
 #env.SharedLibrary('BEM',bem_list)
 
 rough_contact_list = generator_list + surface_list + percolation_list + bem_list
 current_files = set(os.listdir(Dir('.').srcnode().abspath))
 unexpected_files = list(current_files - set(rough_contact_list + ['SConscript']))
 
 unexpected_headers = [ f for f in unexpected_files if os.path.splitext(f)[1] == '.hh' ]
 
 unexpected_files = [ f for f in unexpected_files if not os.path.splitext(f)[1] == '.hh' ]
 
 #print "Unexpected files : {}".format(unexpected_files)
 #print "Include path : {}".format([str(x) for x in env['CPPPATH']])
 #print unexpected_headers
 env.SharedLibrary('Tamaas',rough_contact_list)
diff --git a/src/bem_gigipol.cpp b/src/bem_gigipol.cpp
index f47c8c2..597e12f 100644
--- a/src/bem_gigipol.cpp
+++ b/src/bem_gigipol.cpp
@@ -1,421 +1,422 @@
 /**
  *
  * @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__
 void BemGigipol::computeTractionsFromDisplacements() {
   this->applyInverseInfluenceFunctions(this->true_displacements,this->surface_tractions);
 }
 
 Real BemGigipol::computeEquilibrium(Real epsilon,
                                  Real mean_displacement) {
 
 
     this->computeSpectralInfluenceOverDisplacement();
     this->surface_t = 0.;
     this->surface_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 = SurfaceStatistics::computeAverage(surface, 0);
 
   this->computeGaps();
     this->optimizeToMeanDisplacement(mean_displacement);
 
 std::ofstream file ("output.txt");
   this->computeGaps();
 
 
 
   convergence_iterations.clear();
   nb_iterations = 0;max_iterations=5000;
 
   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);
 
 
       //Projection on admissible space
       this->computeGaps();
     this->enforceMeanDisplacement(mean_displacement);
 
       this->computeGaps();
     this->computePressures();
       f=this->computeStoppingCriterion();
 
       if (nb_iterations % dump_freq ==0){
 
 
         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 << " "<< 'G' << " " << G << " "
                   <<  std::endl;
 
 
 
                 UInt n = surface.size();
                 UInt size = n*n;
                Real crit=0;
                 Real disp_norm=0;
                 for (UInt i = 0; i < size; ++i) {
 
                     crit += this->search_direction(i)*this->search_direction(i);
                          disp_norm += (true_displacements(i)*true_displacements(i));
 
                 }
 
 
                file<< std::scientific << std::setprecision(10)<< nb_iterations << " "
                    <<crit/disp_norm<<  " "<<f<< std::endl;
 
 
 
 
       }
 
       convergence_iterations.push_back(f);
       ++nb_iterations;
 
 
 
 
 
     }
 
     this->computePressures();
 
 
     return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::computeStoppingCriterion() {
 
 
     UInt n = surface.size();
     UInt size = n*n;
 
     Real res = 0;
     Real t_sum = std::abs(SurfaceStatistics::computeSum(this->true_displacements));
     this->computeTractionsFromDisplacements();
     this->functional->computeGradFU();
     const Surface<Real> & gradF = this->functional->getGradF();
     Real min = SurfaceStatistics::computeMinimum(gradF);
 
 
 #pragma omp parallel for reduction(+:res)
     for (UInt i = 0; i < size; ++i) {
         res += std::abs((gradF(i)-min)* (  this->true_displacements(i)-surface(i)));
         // res +=
         //   this->surface_tractions[i].real()
         //   *(surface_displacements[i].real() - surface[i].real());
     }
 
     return res/std::abs(t_sum);
 }
 /* -------------------------------------------------------------------------- */
 
 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() {
   this->computeTractionsFromDisplacements();
   this->functional->computeGradFU();
   const Surface<Real> & gradF = this->functional->getGradF();
   Real min = SurfaceStatistics::computeMinimum(gradF);
   this->surface_tractions -= min;
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
diff --git a/src/surface_generator.cpp b/src/squared_exponential_adhesion_functional.cpp
similarity index 52%
copy from src/surface_generator.cpp
copy to src/squared_exponential_adhesion_functional.cpp
index 287c15e..7f70819 100644
--- a/src/surface_generator.cpp
+++ b/src/squared_exponential_adhesion_functional.cpp
@@ -1,93 +1,93 @@
 /**
  *
- * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
+
+ * @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 <iostream>
-#include <map>
-#include <fstream>
-/* -------------------------------------------------------------------------- */
-#include "surface_generator.hh"
-#include "surface_generator_voss.hh"
-#include "surface_generator_filter.hh"
-#include "surface_generator_filter_bessel.hh"
-#include "surface_generator_filter_fft.hh"
-#include "surface_generator_crenel.hh"
-#include "surface_generator_ellipsoid.hh"
-#include "surface_statistics.hh"
+
+#include "squared_exponential_adhesion_functional.hh"
+#include <cmath>
+
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
-using namespace std;
 /* -------------------------------------------------------------------------- */
-/* Sort functions */
-int sortAscend(const void * a, const void * b) {
-  return ( *(int*)a - *(int*)b );
-}
+
+SquaredExponentialAdhesionFunctional::SquaredExponentialAdhesionFunctional(BemFFTBase & bem):
+  Functional(bem),
+  exponential_term(bem.getSurface().size(), bem.getSurface().getL()),
+  surface_t(bem.getSurface().size(), bem.getSurface().getL())
+{}
+
 /* -------------------------------------------------------------------------- */
 
+SquaredExponentialAdhesionFunctional::~SquaredExponentialAdhesionFunctional()
+{}
 
-int sortDescend(const void * a, const void * b) {
-  return ( *(int*)b - *(int*)a );
-}
 /* -------------------------------------------------------------------------- */
 
-SurfaceGenerator::SurfaceGenerator():
-  rms(0.), random_seed(0),surface(NULL),gridsize(0) {
+Real SquaredExponentialAdhesionFunctional::computeF() {
+  this->computeExponentialTerm();
+  return -SurfaceStatistics::computeSum(this->exponential_term);
 }
 
 /* -------------------------------------------------------------------------- */
 
-void SurfaceGenerator::Init() {
+void SquaredExponentialAdhesionFunctional::computeGradFU() {
+  this->computeExponentialTerm();
 
-  /* Post-process grid size */
-  if (gridsize % 2) {
-    SURFACE_FATAL("Grid size should even and is " << gridsize);
-  }
 
-  surface = new Surface<Real>(gridsize,1.);
 
-  return;
+  *this->gradF += this->exponential_term;
 }
+
 /* -------------------------------------------------------------------------- */
 
-Real SurfaceGenerator::constrainRMS(){
-  int n = surface->size();
-  Real measured_rms = SurfaceStatistics::computeStdev(*surface);
+void SquaredExponentialAdhesionFunctional::computeGradFP() {
+  SURFACE_FATAL("Adhesion functional (exponential) cannot be differentiated with respect to P");
 
-  std::cerr <<  " actual rms is " << rms << endl;
-  /* shift to center to mean value */
-  for (int i = 0 ; i < n ; ++i)
-    for (int j = 0 ; j < n ; ++j){
-      surface->at(i,j) *= rms/measured_rms;
-    }
-  return rms/measured_rms;
 }
 
 /* -------------------------------------------------------------------------- */
 
+void SquaredExponentialAdhesionFunctional::computeExponentialTerm() {
+  const Surface<Real> & gap = bem.getGap();
+  const UInt n = gap.size();
+  const UInt size = n * n;
+
+  Real rho = this->getParameter("rho");
+  Real surface_energy = this->getParameter("surface_energy");
+  if (surface_energy!=0.) {
+#pragma omp parallel for
+
+    for (UInt i = 0 ; i < size ; ++i) {
+      if (gap(i)>=-1.e-8){
+        this->exponential_term(i) = (gap(i))* surface_energy * exp(- (gap(i))*(gap(i)) / 2*rho*rho) / (rho*rho);
+      }
+    }
+  }
+}
 
 __END_TAMAAS__
diff --git a/src/squared_exponential_adhesion_functional.hh b/src/squared_exponential_adhesion_functional.hh
new file mode 100644
index 0000000..d9fe76f
--- /dev/null
+++ b/src/squared_exponential_adhesion_functional.hh
@@ -0,0 +1,68 @@
+/**
+ *
+ * @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/>.
+ *
+ */
+
+/* -------------------------------------------------------------------------- */
+
+#ifndef __SQUARED_EXPONENTIAL_ADHESION_FUNCTIONAL_HH__
+#define __SQUARED_EXPONENTIAL_ADHESION_FUNCTIONAL_HH__
+
+#include "functional.hh"
+#include "bem_fft_base.hh"
+
+__BEGIN_TAMAAS__
+
+/// Energy functional (in displacement) that contains the adhesion term with quadratic gap
+class SquaredExponentialAdhesionFunctional : public Functional {
+public:
+  /// Constructor
+  explicit SquaredExponentialAdhesionFunctional(BemFFTBase & bem);
+
+  /// Destructor
+  virtual ~SquaredExponentialAdhesionFunctional();
+
+public:
+  /// Compute value of objective function
+  virtual Real computeF();
+
+  /// Compute gradient of objective function
+  virtual void computeGradFU();
+
+  /// Compute gradient of objective function
+  virtual void computeGradFP();
+
+  /// Compute the exponential term in functional
+  void computeExponentialTerm();
+
+private:
+  /// Exponential term in functional
+  Surface<Real> exponential_term;
+
+  /// Surface tractions
+  Surface<Real> surface_t;
+};
+
+__END_TAMAAS__
+
+#endif // __SQUARED_EXPONENTIAL_ADHESION_FUNCTIONAL_HH__
diff --git a/src/surface_generator.cpp b/src/surface_generator.cpp
index 287c15e..60806af 100644
--- a/src/surface_generator.cpp
+++ b/src/surface_generator.cpp
@@ -1,93 +1,96 @@
 /**
  *
  * @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 <map>
 #include <fstream>
 /* -------------------------------------------------------------------------- */
 #include "surface_generator.hh"
 #include "surface_generator_voss.hh"
 #include "surface_generator_filter.hh"
 #include "surface_generator_filter_bessel.hh"
 #include "surface_generator_filter_fft.hh"
 #include "surface_generator_crenel.hh"
 #include "surface_generator_ellipsoid.hh"
 #include "surface_statistics.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 using namespace std;
 /* -------------------------------------------------------------------------- */
 /* Sort functions */
 int sortAscend(const void * a, const void * b) {
   return ( *(int*)a - *(int*)b );
 }
 /* -------------------------------------------------------------------------- */
 
 
 int sortDescend(const void * a, const void * b) {
   return ( *(int*)b - *(int*)a );
 }
 /* -------------------------------------------------------------------------- */
 
 SurfaceGenerator::SurfaceGenerator():
   rms(0.), random_seed(0),surface(NULL),gridsize(0) {
 }
 
 /* -------------------------------------------------------------------------- */
 
 void SurfaceGenerator::Init() {
 
   /* Post-process grid size */
   if (gridsize % 2) {
     SURFACE_FATAL("Grid size should even and is " << gridsize);
   }
 
   surface = new Surface<Real>(gridsize,1.);
 
   return;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceGenerator::constrainRMS(){
   int n = surface->size();
   Real measured_rms = SurfaceStatistics::computeStdev(*surface);
 
+  std::cout << "measured rms is "  << measured_rms << std::endl;
+
+
   std::cerr <<  " actual rms is " << rms << endl;
   /* shift to center to mean value */
   for (int i = 0 ; i < n ; ++i)
     for (int j = 0 ; j < n ; ++j){
       surface->at(i,j) *= rms/measured_rms;
     }
   return rms/measured_rms;
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 __END_TAMAAS__
diff --git a/src/surface_generator_filter.cpp b/src/surface_generator_filter.cpp
index 04173a2..78f6537 100644
--- a/src/surface_generator_filter.cpp
+++ b/src/surface_generator_filter.cpp
@@ -1,150 +1,154 @@
 /**
  *
  * @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 "surface_generator_filter.hh"
 #include "fft_plan_manager.hh"
 #include "fftransform.hh"
 #include <cmath>
 #include <iostream>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 SurfaceGeneratorFilter::SurfaceGeneratorFilter(){
   h_coeff = NULL;
   sources = NULL;
   surface_spectral = NULL;
   white_noise_FFT = NULL;
   h_coeff_power = NULL;
 }
 
 /* -------------------------------------------------------------------------- */
 
 SurfaceGeneratorFilter::~SurfaceGeneratorFilter(){
   delete white_noise_FFT;
   delete surface_spectral;
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 void SurfaceGeneratorFilter::Init(){
   SurfaceGenerator::Init();
 }
 
 /* -------------------------------------------------------------------------- */
 
 void SurfaceGeneratorFilter::applyFilterOnSource(){
   std::cerr << "resize surface" << std::endl;
   if (surface_spectral == NULL)
     surface_spectral = new SurfaceComplex<Real>(surface->size(),surface->getL());
   else surface_spectral->setGridSize(surface->size());
   std::cerr << "do full product in fourier space" << std::endl;
   *surface_spectral = *white_noise_FFT;
   *surface_spectral *= *h_coeff_power;
   std::cerr << "fourier transform filter spectrum" << std::endl;
 
   FFTransform<Real, 2> & plan = FFTPlanManager::get().createPlan(*surface, *surface_spectral);
   plan.backwardTransform();
   FFTPlanManager::get().destroyPlan(*surface, *surface_spectral);
 }
 /* -------------------------------------------------------------------------- */
 
 
 void SurfaceGeneratorFilter::generateWhiteNoiseFFT(){
   white_noise_FFT = new SurfaceComplex<Real>(sources->size(), sources->getL());
   FFTransform<Real, 2> & plan = FFTPlanManager::get().createPlan(*sources, *white_noise_FFT);
   plan.forwardTransform();
   FFTPlanManager::get().destroyPlan(*sources, *white_noise_FFT);
   // white_noise_FFT->dumpToParaview("white noise FFT",NO_REPLICATION);
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 void SurfaceGeneratorFilter::generateHFFT(){
   FFTransform<Real, 2> & plan = FFTPlanManager::get().createPlan(*h_coeff, *h_coeff_power);
   // fourier transform it
   plan.forwardTransform();
   // make it real so that applicable filter is constructed
   h_coeff_power->makeItRealByAbs();
   FFTPlanManager::get().destroyPlan(*h_coeff, *h_coeff_power);
 }
 /* -------------------------------------------------------------------------- */
 
 
 Surface<Real> & SurfaceGeneratorFilter::buildSurface(){
 
   if (surface == NULL) SURFACE_FATAL("Init function was not called");
 
   int n = surface->size();
 
   if (white_noise_FFT == NULL){
     // create the random source
     sources = new Surface<Real>(n,1.);
     Surface<Real>::generateWhiteNoiseSurface(*sources,random_seed);
     fprintf(stderr,"build FFT of white noise\n");
     generateWhiteNoiseFFT();
     delete (sources);
   }
 
   fprintf(stderr,"build filter coefficients\n");
   computeFilterCoefficients();
   if (h_coeff_power == NULL) {
     h_coeff_power = new SurfaceComplex<Real>(h_coeff->size(), h_coeff->getL());
     generateHFFT();
   }
   //  s->setHeights(h_coeff_power);
   //  h_coeff_power->dumpToParaview("h_coeff_power",NO_REPLICATION);
 
   fprintf(stderr,"apply filter on white noise\n");
   applyFilterOnSource();
 
   delete h_coeff_power;
 
   fprintf(stderr,"recenter arround mean the surface\n");
   surface->recenterHeights();
+
+
+
+
   constrainRMS();
 
   //  s->setHeights(heights);
   //s->SetHeights(sources);
   //s->SetHeights(h);
 
 
   return *surface;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void SurfaceGeneratorFilter::clearSource(){
   if (white_noise_FFT != NULL){
     delete white_noise_FFT;
     white_noise_FFT = NULL;
   }
 }
 
 __END_TAMAAS__