diff --git a/python/bem.i b/python/bem.i
index 44fbea5..3e654f7 100644
--- a/python/bem.i
+++ b/python/bem.i
@@ -1,50 +1,51 @@
 /**
  *
  * @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 "functional.hh"
 %include "meta_functional.hh"
 %include "elastic_energy_functional.hh"
 %include "complimentary_term_functional.hh"
 %include "exponential_adhesion_functional.hh"
+%include "maugis_adhesion_functional.hh"
 %include "fftransform.hh"
 
 namespace tamaas {
 %template(FFTransformReal) FFTransform<Real>;
 }
diff --git a/python/surface.i b/python/surface.i
index 02650c9..8eb084d 100644
--- a/python/surface.i
+++ b/python/surface.i
@@ -1,257 +1,258 @@
 /**
  *
  * @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"
 #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 "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 "maugis_adhesion_functional.hh"
 #include "surface_timer.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>)
 
  /* -------------------------------------------------------------------------- */
 /* Wrapping Surface class                                                     */
 /* -------------------------------------------------------------------------- */
 
 %include "tamaas.hh"
 
 %ignore tamaas::Grid::operator+=;
 %ignore tamaas::Grid::operator*=;
 %ignore tamaas::Grid::operator-=;
 %ignore tamaas::Grid::operator/=;
 
 %include "grid.hh"
 namespace tamaas {
   %template (Grid2dReal) Grid<Real, 2>;
   %template (Grid2dComplex) Grid<complex, 2>;
   %template (Grid2dInt) Grid<Int, 2>;
 }
 
 %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 <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 %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{
     UInt _n = PyArray_NDIM((PyArrayObject*)$input);
     if (_n != 2) SURFACE_FATAL("incompatible numpy dimension " << _n);
     npy_intp * ndims = PyArray_DIMS((PyArrayObject*)$input);
     UInt sz = ndims[0];
     UInt sz2 = ndims[1];
     if (sz != sz2) SURFACE_FATAL("incompatible numpy dimensions " << sz << " " << sz2);
     PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew($input);
     if (iter == NULL) throw;
     $1 = new SurfaceForPython< DATA_TYPE >(( DATA_TYPE *)(iter->dataptr),sz,1.);
     PyArray_ITER_NEXT(iter);
   }
  }
 
 %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  )
 
 namespace tamaas {
 %template(SurfaceForPythonReal) SurfaceForPython< Real >;
 }
 
 %include "surface_timer.hh"
diff --git a/src/SConscript b/src/SConscript
index ead2494..1048949 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,71 +1,72 @@
 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
+maugis_adhesion_functional.cpp
 complimentary_term_functional.cpp
 fftransform.cpp
 fftransform_fftw.cpp
 fft_plan_manager.cpp
 grid.cpp
 types.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 18692dd..04d218d 100644
--- a/src/bem_gigipol.cpp
+++ b/src/bem_gigipol.cpp
@@ -1,396 +1,396 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #include <vector>
 #include "surface.hh"
 #include "bem_gigipol.hh"
 #include <iostream>
 #include <sstream>
 #include <fstream>
 #include <iomanip>
 #include <sstream>
 #include <cmath>
 /* -------------------------------------------------------------------------- */
 #define TIMER
 #include "surface_timer.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 Real BemGigipol::computeEquilibrium(Real epsilon,
                                  Real mean_displacement) {
 
 
     this->computeSpectralInfluenceOverDisplacement();
     this->surface_t = 0.;
     this->surface_r = 0.;
     this->surface_spectral_r = 0.;
     this->search_direction = 0.;
     this->pold = 0.;
     this->surface_displacements = 0.;
     Real delta = 0.;
     Real Gold = 1.;
     Real f = 1e300;
     Real fPrevious = 1e300;
     UInt n = surface.size();
     UInt size = n*n;
 
-// this->true_displacements =0;
+ this->true_displacements =0;
  //this->true_displacements = SurfaceStatistics::computeAverage(surface, 0);
 
   this->computeGaps();
     this->optimizeToMeanDisplacement(mean_displacement);
 
 
   this->computeGaps();
 
   convergence_iterations.clear();
-  nb_iterations = 0;max_iterations=10000;
+  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;
 
       }
 
       convergence_iterations.push_back(f);
       ++nb_iterations;
 
 
 
 
 
     }
 
     this->computePressures();
 
 
     return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigipol::computeStoppingCriterion() {
 
 
-    Real rho = this->functional->getParameter("rho");
-    Real surface_energy = this->functional->getParameter("surface_energy");
-
-     Real crit=0.;
-
-Real norm = std::abs(SurfaceStatistics::computeSum(this->true_displacements));
-     UInt n = surface.size();
-     UInt size = n*n;
+    UInt n = surface.size();
+    UInt size = n*n;
 
-   #pragma omp parallel for reduction(+:crit)
-     for (UInt i = 0; i < size; ++i) {
-      crit += std::abs(this->gap(i)*(this->surface_tractions(i)+surface_energy/rho * exp(- (gap(i)) / rho)));
+    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);
 
-//crit += 1;
-     }
 
+#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 crit/norm ;
+    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/bem_kato.cpp b/src/bem_kato.cpp
index 994d137..1d00b90 100644
--- a/src/bem_kato.cpp
+++ b/src/bem_kato.cpp
@@ -1,301 +1,300 @@
 /**
  *
  * @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_kato.hh"
 #include <iostream>
 #include <sstream>
 #include <fstream>
 #include <iomanip>
 #include <sstream>
 #include <cmath>
 /* -------------------------------------------------------------------------- */
 #define TIMER
 #include "surface_timer.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 
 Real BemKato::linescan(Real scale,Real pressure){
   updatePressure(scale);
   shiftPressure(pressure);
   truncatePressure();
   Real res = computeF();
 
   return res;
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real BemKato::linesearch(Real & hmax, Real fold,
 			Real pressure, int search_flag){
   if (search_flag){
     Real h  = hmax;
     //  Real fold = bem.computeF();
     //if (fold == 0) fold = 1e300;
     Real f = linescan(h,pressure);
     if (f < fold) return f;
 
     while (f > fold){
       h *= 0.5;
       if (h < 1e-3) throw 1;
       f = linescan(h,pressure);
     }
     f = linescan(h,pressure);
     //  if (hmax / h > 10) hmax /=2;
     return f;
   }
   return linescan(hmax,pressure);
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 Real BemKato::computeEquilibrium(Real epsilon,
 				Real pressure){
 
   //  UInt n = surface.size();
   //  UInt size = n*n;
 
   this->computeSpectralInfluenceOverDisplacement();
   this->computeDisplacementsFromTractions();
   this->functional->computeGradFP();
   this->backupTractions();
   Real f = 1.;
   Real fPrevious = 1e300;
 
   this->nb_iterations = 0;
   this->convergence_iterations.clear();
 
   while(f > epsilon && this->nb_iterations < this->max_iterations) {
 
     fPrevious = f;
     this->computeDisplacementsFromTractions();
     this->functional->computeGradFP();
     this->backupTractions();
     try {
       f = linescan(1.,pressure);
     }
     catch (int e){
       std::cout << " line search failed " << std::endl;
       f = linescan(1,pressure);
       nb_iterations = 0;
       break;
     }
     if (nb_iterations % dump_freq == 0) {
       //    std::cout << std::scientific << std::setprecision(10) << nb_iterations << " " << f << " " << f-fold << " " << ((f-fold)/forigin) <<  std::endl;
       std::cout << std::scientific << std::setprecision(10)
                 << nb_iterations << " "
                 << f << " " << f-fPrevious
                 << std::endl;
     }
 
     convergence_iterations.push_back(f);
     ++nb_iterations;
   }
 
 
   this->computeTrueDisplacements();
   this->computeGaps();
 
 
 
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 void BemKato::updatePressure(Real scale){
   STARTTIMER("updatePressure");
   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->surface_tractions(i) =
       this->surface_tractions_backup(i) - gradF(i) * scale;
   }
   STOPTIMER("updatePressure");
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemKato::backupTractions(){
   STARTTIMER("switchPressure");
 
   UInt n = surface.size();
   UInt size = n*n;
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     this->surface_tractions_backup(i) = this->surface_tractions(i);
   }
 
   STOPTIMER("switchPressure");
 }
 
 
 /* -------------------------------------------------------------------------- */
 Real BemKato::positivePressure(Real step){
   STARTTIMER("positivePressure");
 
   UInt n = surface.size();
   UInt size = n*n;
   Real p_tot = 0.0;
 
 #pragma omp parallel for reduction(+:p_tot)
   for (UInt i = 0; i < size; ++i) {
     Real sh_press = this->surface_tractions(i) + step;
     if (sh_press > max_pressure) p_tot += max_pressure;
     else p_tot += sh_press*(sh_press > 0);
   }
 
   STOPTIMER("positivePressure");
   return p_tot/n/n;
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 void BemKato::shiftPressure(Real required_pressure){
 
   Real step_min = -10;
   Real step_max =  10;
 
   STARTTIMER("shiftPressureInitialGuess");
 
   Real p_max = positivePressure(step_max);
   Real p_min = positivePressure(step_min);
 
 
   for(UInt i = 0; p_max < required_pressure && i < 8; ++i, step_max*=10) {
     p_max = positivePressure(step_max);
   }
   for(UInt i = 0; p_min > required_pressure && i < 8; ++i, step_min*=10) {
     p_min = positivePressure(step_min);
   }
 
   Real p = positivePressure(0.0);
   Real epsilon = 1e-12;
 
   STOPTIMER("shiftPressureInitialGuess");
 
   STARTTIMER("shiftPressureDichotomy");
   while (fabs(step_min - step_max) > epsilon){
     Real step = (step_min + step_max)/2;
     p = positivePressure(step);
     if (p > required_pressure) step_max = step;
     else if (p < required_pressure) step_min = step;
     else {
       step_max = step;
       step_min = step;
     }
   }
   STOPTIMER("shiftPressureDichotomy");
 
   STARTTIMER("shiftPressureTrueShift");
   UInt n = surface.size();
   UInt size = n*n;
   //shift the pressure so that satisfies the constraint
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     this->surface_tractions(i) += (step_max+step_min)/2;
   }
   STOPTIMER("shiftPressureTrueShift");
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemKato::truncatePressure(){
   STARTTIMER("truncatePressure");
 
   UInt n = surface.size();
   UInt size = n*n;
   //shift the pressure so that satisfies the constraint
 #pragma omp parallel for
   for (UInt i = 0; i < size; ++i) {
     if (this->surface_tractions(i) > max_pressure)
       this->surface_tractions(i) = max_pressure;
 
     else this->surface_tractions(i) *= (this->surface_tractions(i) > 0);
   }
 
   STOPTIMER("truncatePressure");
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 
 void BemKato::computeTrueDisplacements() {
 
   UInt n = surface.size();
   UInt size = n*n;
 
   Real shift = 1e300;
 
 #pragma omp parallel for reduction(min: shift)
   for (UInt i = 0 ; i < size ; i++) {
     if (surface_displacements(i) - shift - surface(i) < 0. &&
 	surface_tractions(i) < max_pressure){
       shift = surface_displacements(i) - surface(i);
     }
   }
 
   this->true_displacements = surface_displacements;
   this->true_displacements -= shift;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemKato::computeF() {
   UInt n = surface.size();
   UInt size = n*n;
 
   Real res = 0;
-  Real d_max = SurfaceStatistics::computeMaximum(surface_displacements);
-  Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions));
-  Real s_max = SurfaceStatistics::computeMaximum(surface);
 
+  Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions));
+  computeTrueDisplacements();
 #pragma omp parallel for reduction(+:res)
   for (UInt i = 0; i < size; ++i) {
     if (this->surface_tractions(i) == this->max_pressure) continue;
-    res += std::abs(surface_tractions(i) * (surface_displacements(i)-d_max-surface(i)+s_max));
+    res += std::abs(surface_tractions(i) * (  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);
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 __END_TAMAAS__
diff --git a/src/bem_polonski.cpp b/src/bem_polonski.cpp
index 17f23c7..64e4c2e 100644
--- a/src/bem_polonski.cpp
+++ b/src/bem_polonski.cpp
@@ -1,437 +1,436 @@
 /**
  *
  * @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_polonski.hh"
 #include "fft_plan_manager.hh"
 #include "fftransform.hh"
 #include <iostream>
 #include <sstream>
 #include <fstream>
 #include <iomanip>
 #include <sstream>
 #include <cmath>
 /* -------------------------------------------------------------------------- */
 #define TIMER
 #include "surface_timer.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 BemPolonski::BemPolonski(Surface<Real> & p) :
-  BemFFTBase(p),
-  surface_t(p.size(),p.getL()),
-  surface_r(p.size(),p.getL()),
-  pold(p.size(),p.getL()),
-  p_sat(-1.){
-  e_star = 1.;
-  max_iterations = 2000;
+    BemFFTBase(p),
+    surface_t(p.size(),p.getL()),
+    surface_r(p.size(),p.getL()),
+    pold(p.size(),p.getL()),
+    p_sat(-1.){
+    e_star = 1.;
+    max_iterations = 2000;
 }
 
 /* -------------------------------------------------------------------------- */
 BemPolonski::~BemPolonski() {}
 /* -------------------------------------------------------------------------- */
 
 Real BemPolonski::computeEquilibrium(Real epsilon,
                                      Real pressure) {
-  this->computeSpectralInfluenceOverDisplacement();
-  this->surface_t = 0.;
-  this->surface_r = 0.;
-  this->pold = 0.;
-  this->surface_displacements = 0.;
-  Real delta = 0.;
-  Real Gold = 1.;
-  Real f = 1e300;
-  Real fPrevious = 1e300;
+    this->computeSpectralInfluenceOverDisplacement();
+    this->surface_t = 0.;
+    this->surface_r = 0.;
+    this->pold = 0.;
+    this->surface_displacements = 0.;
+    Real delta = 0.;
+    Real Gold = 1.;
+    Real f = 1e300;
+    Real fPrevious = 1e300;
+
+    Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions);
+    if (current_pressure <= 0.) surface_tractions = pressure;
+    this->enforcePressureBalance(pressure);
+    convergence_iterations.clear();
+    nb_iterations = 0;
 
-  Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions);
-  if (current_pressure <= 0.) surface_tractions = pressure;
-  this->enforcePressureBalance(pressure);
-  convergence_iterations.clear();
-  nb_iterations = 0;
+    while(f > epsilon && nb_iterations < max_iterations) {
+        fPrevious = f;
+        this->computeDisplacementsFromTractions();
+        this->computeGaps();
 
-  while(f > epsilon && nb_iterations < max_iterations) {
-    fPrevious = f;
-    this->computeDisplacementsFromTractions();
-    this->computeGaps();
+        Real gbar = this->computeMeanGapsInContact();
 
-    Real gbar = this->computeMeanGapsInContact();
+        this->gap -= gbar;
+        Real G = this->computeG();
+        std::cout << std::scientific << std::setprecision(10)    << G<<  std::fixed << std::endl;
 
-    this->gap -= gbar;
-    Real G = this->computeG();
-    // std::cout << std::scientific << std::setprecision(10)    << G<<  std::fixed << std::endl;
 
+        this->updateT(G,Gold,delta);
+        Real tau = this->computeTau();
+        pold = this->surface_tractions;
+        Gold = G;
+        delta = this->updateTractions(tau);
 
-    this->updateT(G,Gold,delta);
-    Real tau = this->computeTau();
-    pold = this->surface_tractions;
-    Gold = G;
-    delta = this->updateTractions(tau);
 
+        //Projection on admissible space
 
-    //Projection on admissible space
+        this->enforcePressureBalance(pressure);
+        // std::cout << std::scientific << std::setprecision(10)
+        //     << SurfaceStatistics::computeAverage(surface_tractions) <<  std::fixed << std::endl;
+        //
 
-    this->enforcePressureBalance(pressure);
-   // std::cout << std::scientific << std::setprecision(10)
-         //     << SurfaceStatistics::computeAverage(surface_tractions) <<  std::fixed << std::endl;
-//
+        f = this->computeF();
 
-    f = this->computeF();
+        if (nb_iterations % dump_freq == 0){
+            Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions);
 
-    if (nb_iterations % dump_freq == 0){
-      Real A = SurfaceStatistics::computeContactAreaRatio(surface_tractions);
+            std::cout << std::scientific << std::setprecision(10)
+                      << nb_iterations << " "
+                      << f << " " << f-fPrevious << " " << A
+                      << std::fixed << std::endl;
+        }
 
-      std::cout << std::scientific << std::setprecision(10)
-                << nb_iterations << " "
-                << f << " " << f-fPrevious << " " << A
-                << std::fixed << std::endl;
+        convergence_iterations.push_back(f);
+        ++nb_iterations;
     }
-
-    convergence_iterations.push_back(f);
-    ++nb_iterations;
-  }
-  this->computeTrueDisplacements();
-  return f;
+    this->computeTrueDisplacements();
+    return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemPolonski::computeGaps() {
-  UInt size = surface.size();
+    UInt size = surface.size();
 
 #pragma omp parallel for
-  for (UInt i=0; i < size*size; i++) {
-    gap(i) = surface_displacements(i) - surface(i);
-  }
+    for (UInt i=0; i < size*size; i++) {
+        gap(i) = surface_displacements(i) - surface(i);
+    }
 }
 
 
 /* -------------------------------------------------------------------------- */
 Real BemPolonski::computeMeanGapsInContact(){
-  STARTTIMER("computeMeanGapsInContact");
+    STARTTIMER("computeMeanGapsInContact");
 
-  UInt n = surface.size();
-  UInt size = n*n;
-  Real res = 0.;
-  UInt nb_contact = 0;
+    UInt n = surface.size();
+    UInt size = n*n;
+    Real res = 0.;
+    UInt nb_contact = 0;
 #pragma omp parallel for reduction(+: nb_contact,res)
-  for (UInt i = 0; i < size; ++i) {
-    if (this->surface_tractions(i) <= 0 ||std::abs(this->surface_tractions(i) - p_sat)< 1.0e-10 ) continue;
-    ++nb_contact;
-    res += this->gap(i);
-  }
-
-  res /= nb_contact;
-  STOPTIMER("computeMeanGapsInContact");
-  return res;
+    for (UInt i = 0; i < size; ++i) {
+        if (this->surface_tractions(i) <= 0 ||std::abs(this->surface_tractions(i) - p_sat)< 1.0e-10 ) continue;
+        ++nb_contact;
+        res += this->gap(i);
+    }
+
+    res /= nb_contact;
+    STOPTIMER("computeMeanGapsInContact");
+    return res;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemPolonski::computeG(){
-  STARTTIMER("computeG");
+    STARTTIMER("computeG");
 
-  UInt n = surface.size();
-  UInt size = n*n;
-  Real res = 0.;
+    UInt n = surface.size();
+    UInt size = n*n;
+    Real res = 0.;
 #pragma omp parallel for reduction(+: res)
-  for (UInt i = 0; i < size; ++i) {
-    if (this->surface_tractions(i) <= 0   ||std::abs(this->surface_tractions(i) - p_sat)< 1.0e-10  ) continue;
-    Real val = this->gap(i);
-    res += val*val;
-  }
-
-  STOPTIMER("computeG");
-  return res;
+    for (UInt i = 0; i < size; ++i) {
+        if (this->surface_tractions(i) <= 0   ||std::abs(this->surface_tractions(i) - p_sat)< 1.0e-10  ) continue;
+        Real val = this->gap(i);
+        res += val*val;
+    }
+
+    STOPTIMER("computeG");
+    return res;
 }
 /* -------------------------------------------------------------------------- */
 
 
 void BemPolonski::updateT(Real G, Real Gold, Real delta){
-  STARTTIMER("updateT");
+    STARTTIMER("updateT");
 
-  UInt n = surface.size();
-  UInt size = n*n;
-  Real factor = delta*G/Gold;
+    UInt n = surface.size();
+    UInt size = n*n;
+    Real factor = delta*G/Gold;
 #pragma omp parallel for
-  for (UInt i = 0; i < size; ++i) {
-    if (this->surface_tractions(i) <= 0) this->surface_t(i) = 0.;
-    else{
-      this->surface_t(i) *= factor;
-      this->surface_t(i) += this->gap(i);
+    for (UInt i = 0; i < size; ++i) {
+        if (this->surface_tractions(i) <= 0) this->surface_t(i) = 0.;
+        else{
+            this->surface_t(i) *= factor;
+            this->surface_t(i) += this->gap(i);
+        }
     }
-  }
 
-  STOPTIMER("updateT");
+    STOPTIMER("updateT");
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 Real BemPolonski::computeTau(){
-  STARTTIMER("computeTau");
+    STARTTIMER("computeTau");
 
-  this->applyInfluenceFunctions(surface_t, surface_r);
+    this->applyInfluenceFunctions(surface_t, surface_r);
 
-  const UInt size = surface_t.size() * surface_t.size();
+    const UInt size = surface_t.size() * surface_t.size();
 
-  Real rbar = 0;
-  UInt nb_contact = 0;
+    Real rbar = 0;
+    UInt nb_contact = 0;
 #pragma omp parallel for reduction(+: nb_contact,rbar)
-  for (UInt i = 0; i < size; ++i) {
-    if (this->surface_tractions(i) <= 0) continue;
-    ++nb_contact;
-    rbar += surface_r(i);
-  }
-  rbar /= nb_contact;
-  surface_r -= rbar;
+    for (UInt i = 0; i < size; ++i) {
+        if (this->surface_tractions(i) <= 0) continue;
+        ++nb_contact;
+        rbar += surface_r(i);
+    }
+    rbar /= nb_contact;
+    surface_r -= rbar;
 
-  Real tau_sum1 = 0.,tau_sum2 = 0.;
+    Real tau_sum1 = 0.,tau_sum2 = 0.;
 
 #pragma omp parallel for reduction(+: tau_sum1, tau_sum2)
-  for (UInt i = 0; i < size; ++i) {
-    if (this->surface_tractions(i) <= 0) continue;
-    tau_sum1 += this->gap(i)*surface_t(i);
-    tau_sum2 += surface_r(i)*surface_t(i);
-  }
-  Real tau = tau_sum1/tau_sum2;
-
-  STOPTIMER("computeTau");
-  return tau;
+    for (UInt i = 0; i < size; ++i) {
+        if (this->surface_tractions(i) <= 0) continue;
+        tau_sum1 += this->gap(i)*surface_t(i);
+        tau_sum2 += surface_r(i)*surface_t(i);
+    }
+    Real tau = tau_sum1/tau_sum2;
+
+    STOPTIMER("computeTau");
+    return tau;
 }
 
 
 
 /* -------------------------------------------------------------------------- */
 
 Real BemPolonski::updateTractions(Real tau){
-   STARTTIMER("updateTractions");
+    STARTTIMER("updateTractions");
 
-   UInt n = surface.size();
-   UInt size = n*n;
+    UInt n = surface.size();
+    UInt size = n*n;
 #pragma omp parallel for
-   for (UInt i = 0; i < size; ++i) {
-     if (this->surface_tractions(i) <= 0) continue;
-     this->surface_tractions(i) -= tau*this->surface_t(i);
-     if (this->surface_tractions(i) < 0.) this->surface_tractions(i) = 0;
-   }
-
-   //compute number of interpenetration without contact
-   UInt nb_iol = 0;
+    for (UInt i = 0; i < size; ++i) {
+        if (this->surface_tractions(i) <= 0) continue;
+        this->surface_tractions(i) -= tau*this->surface_t(i);
+        if (this->surface_tractions(i) < 0.) this->surface_tractions(i) = 0;
+    }
+
+    //compute number of interpenetration without contact
+    UInt nb_iol = 0;
 #pragma omp parallel for reduction(+: nb_iol)
-   for (UInt i = 0; i < size; ++i) {
-     if (this->surface_tractions(i) <= 0. &&
-         this->gap(i) < 0.) {
-       this->surface_tractions(i) -= tau*this->gap(i);
-       ++nb_iol;
-     }
-   }
+    for (UInt i = 0; i < size; ++i) {
+        if (this->surface_tractions(i) <= 0. &&
+                this->gap(i) < 0.) {
+            this->surface_tractions(i) -= tau*this->gap(i);
+            ++nb_iol;
+        }
+    }
 
 
-   //compute number of interpenetration without contact
-   UInt nb_sl = 0;
+    //compute number of interpenetration without contact
+    UInt nb_sl = 0;
 #pragma omp parallel for reduction(+: nb_sl)
-   for (UInt i = 0; i < size; ++i) {
-     if (this->surface_tractions(i) ==p_sat &&
-        0.<  this->gap(i) ) {
-       this->surface_tractions(i) -= tau*this->gap(i);
-       ++nb_sl;
-     }
-   }
+    for (UInt i = 0; i < size; ++i) {
+        if (this->surface_tractions(i) ==p_sat &&
+                0.<  this->gap(i) ) {
+            this->surface_tractions(i) -= tau*this->gap(i);
+            ++nb_sl;
+        }
+    }
 
 
 
 
-   Real delta = 0;
-   if (nb_iol > 0 || nb_sl> 0) delta = 0.;
-   else            delta = 1.;
+    Real delta = 0;
+    if (nb_iol > 0 || nb_sl> 0) delta = 0.;
+    else            delta = 1.;
 
-   STOPTIMER("updateTractions");
-   return delta;
+    STOPTIMER("updateTractions");
+    return delta;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemPolonski::enforcePressureBalance(Real applied_pressure){
-   STARTTIMER("enforcePressureBalance");
+    STARTTIMER("enforcePressureBalance");
 
-   UInt n = surface.size();
-   UInt size = n*n;
+    UInt n = surface.size();
+    UInt size = n*n;
 
 
-   // If we have no saturation, scale to match applied_pressure
-   if (this->p_sat <= 0. ||
-       SurfaceStatistics::computeMaximum(this->surface_tractions) < this->p_sat) {
-    Real pressure = 0.;
+    // If we have no saturation, scale to match applied_pressure
+    if (this->p_sat <= 0. ||
+            SurfaceStatistics::computeMaximum(this->surface_tractions) < this->p_sat) {
+        Real pressure = 0.;
 #pragma omp parallel for reduction(+: pressure)
-     for (UInt i = 0; i < size; ++i) {
-      pressure += this->surface_tractions(i);
-    }
-
-     pressure *= 1./size;
-    this->surface_tractions *= applied_pressure/pressure;
-  }
-
-  // If we have saturation, use secant method to find zero of saturationFunction
-   else {
-
-
-
-
-//   // Checking if we can find a zero
-//    Real limit =
-//      this->p_sat * SurfaceStatistics::computeContactArea(this->surface_tractions)
-//       - applied_pressure;
-//     if (limit < 0) SURFACE_FATAL("Saturation pressure is too small for applied pressure");
-
-//     // Initial points
-//     Real x_n_2 = 0., x_n_1 = 1., x_n = 0.;
-//     Real f_n_2 = 0., f_n_1 = 0., f_n = 0.;
-
-//     // Secant loop
-//     do {
-//       f_n_2 = saturationFunction(x_n_2, applied_pressure);
-//       f_n_1 = saturationFunction(x_n_1, applied_pressure);
-//       x_n = x_n_1 - f_n_1 * (x_n_1 - x_n_2) / (f_n_1 - f_n_2);
-//       f_n = saturationFunction(x_n, applied_pressure);
-//      x_n_2 = x_n_1;
-//       x_n_1 = x_n;
-//    } while (std::abs(f_n) > 1e-16);
-
-//    this->surface_tractions *= x_n;
-//     // Truncating pressures
-//#pragma omp parallel for
-//     for (UInt i = 0 ; i < size ; i++) {
-//      if (this->surface_tractions(i) > this->p_sat)
-//         this->surface_tractions(i) = this->p_sat;
-//     }
-
-     Real mu_I_sat=0.;
-             for (UInt i = 0; i < size; ++i) {
-                 if ( this->p_sat< this->surface_tractions(i)   ){
-              mu_I_sat+=1.;
-                 }
-             }
+        for (UInt i = 0; i < size; ++i) {
+            pressure += this->surface_tractions(i);
+        }
 
+        pressure *= 1./size;
+        this->surface_tractions *= applied_pressure/pressure;
+    }
 
-             Real pc=0.;
+    // If we have saturation, use secant method to find zero of saturationFunction
+    else {
+
+
+
+
+        //   // Checking if we can find a zero
+        //    Real limit =
+        //      this->p_sat * SurfaceStatistics::computeContactArea(this->surface_tractions)
+        //       - applied_pressure;
+        //     if (limit < 0) SURFACE_FATAL("Saturation pressure is too small for applied pressure");
+
+        //     // Initial points
+        //     Real x_n_2 = 0., x_n_1 = 1., x_n = 0.;
+        //     Real f_n_2 = 0., f_n_1 = 0., f_n = 0.;
+
+        //     // Secant loop
+        //     do {
+        //       f_n_2 = saturationFunction(x_n_2, applied_pressure);
+        //       f_n_1 = saturationFunction(x_n_1, applied_pressure);
+        //       x_n = x_n_1 - f_n_1 * (x_n_1 - x_n_2) / (f_n_1 - f_n_2);
+        //       f_n = saturationFunction(x_n, applied_pressure);
+        //      x_n_2 = x_n_1;
+        //       x_n_1 = x_n;
+        //    } while (std::abs(f_n) > 1e-16);
+
+        //    this->surface_tractions *= x_n;
+        //     // Truncating pressures
+        //#pragma omp parallel for
+        //     for (UInt i = 0 ; i < size ; i++) {
+        //      if (this->surface_tractions(i) > this->p_sat)
+        //         this->surface_tractions(i) = this->p_sat;
+        //     }
+        Real mu_I_sat=0.;
+#pragma omp parallel for reduction(+: mu_I_sat)
+        for (UInt i = 0; i < size; ++i) {
+            if ( this->p_sat< this->surface_tractions(i)   ){
+                mu_I_sat+=1.;
+            }
+        }
+
+
+        Real pc=0.;
 #pragma omp parallel for reduction(+: pc)
- for (UInt i = 0; i < size; ++i) {
-     if ( this->surface_tractions(i) < this->p_sat  ){
- pc+=this->surface_tractions(i);}
- }
+        for (UInt i = 0; i < size; ++i) {
+            if ( this->surface_tractions(i) < this->p_sat  ){
+                pc+=this->surface_tractions(i);}
+        }
 
 #pragma omp parallel for
- for (UInt i = 0; i < size; ++i) {
-     if ( this->surface_tractions(i) < this->p_sat   ){
- this->surface_tractions(i)= this->surface_tractions(i)*((size*applied_pressure-this->p_sat*mu_I_sat)/pc);}
- }
+        for (UInt i = 0; i < size; ++i) {
+            if ( this->surface_tractions(i) < this->p_sat   ){
+                this->surface_tractions(i)= this->surface_tractions(i)*((size*applied_pressure-this->p_sat*mu_I_sat)/pc);}
+        }
 #pragma omp parallel for
- for (UInt i = 0; i < size; ++i) {
-     if ( this->p_sat< this->surface_tractions(i)   ){
-this->surface_tractions(i)=this->p_sat;
-     }
- }
+        for (UInt i = 0; i < size; ++i) {
+            if ( this->p_sat< this->surface_tractions(i)   ){
+                this->surface_tractions(i)=this->p_sat;
+            }
+        }
 
 
-  }
+    }
 
 
 
 
 
-   STOPTIMER("enforcePressureBalance");
- }
+    STOPTIMER("enforcePressureBalance");
+}
 
 /* -------------------------------------------------------------------------- */
 
 Real BemPolonski::saturationFunction(Real alpha, Real applied_pressure) {
-  const UInt n = surface.size();
-  const UInt size = n*n;
-  Real sum = 0.;
+    const UInt n = surface.size();
+    const UInt size = n*n;
+    Real sum = 0.;
 #pragma omp parallel for reduction(+: sum)
-  for (UInt i = 0 ; i < size ; i++) {
-    Real alpha_p = alpha * this->surface_tractions(i);
-    Real alpha_p_I = (alpha_p > this->p_sat)? alpha_p - this->p_sat : 0.;
-    sum += alpha_p - alpha_p_I;
-  }
-  sum /= size;
-  return sum - applied_pressure;
+    for (UInt i = 0 ; i < size ; i++) {
+        Real alpha_p = alpha * this->surface_tractions(i);
+        Real alpha_p_I = (alpha_p > this->p_sat)? alpha_p - this->p_sat : 0.;
+        sum += alpha_p - alpha_p_I;
+    }
+    sum /= size;
+    return sum - applied_pressure;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemPolonski::computeF() {
-  UInt n = surface.size();
-  UInt size = n*n;
+    UInt n = surface.size();
+    UInt size = n*n;
 
-  Real res = 0;
-  Real d_max = SurfaceStatistics::computeMaximum(surface_displacements);
-  Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions));
-  Real s_max = SurfaceStatistics::computeMaximum(surface);
+    Real res = 0;
 
+    Real t_sum = std::abs(SurfaceStatistics::computeSum(surface_tractions));
+    computeTrueDisplacements();
 #pragma omp parallel for reduction(+:res)
-  for (UInt i = 0; i < size; ++i) {
-    res += std::abs(surface_tractions(i) * (surface_displacements(i)-d_max-surface(i)+s_max));
-    // res +=
-    //   this->surface_tractions[i].real()
-    //   *(surface_displacements[i].real() - surface[i].real());
-  }
-
-  return res/t_sum;
+    for (UInt i = 0; i < size; ++i) {
+        if (this->surface_tractions(i) == this->p_sat) continue;
+        res += std::abs(surface_tractions(i) * (  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);
 }
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 
 void BemPolonski::computeTrueDisplacements(){
     this->applyInfluenceFunctions(this->surface_tractions, this->surface_displacements);
 
-  this->true_displacements = this->surface_displacements;
+    this->true_displacements = this->surface_displacements;
 
-  UInt n = surface.size();
-  UInt size = n*n;
+    UInt n = surface.size();
+    UInt size = n*n;
 
-  Real shift = 1e300;
+    Real shift = 1e300;
 
 
+#pragma omp parallel for reduction(min:shift)
+    for (UInt i = 0; i < size; ++i) {
 
-  for (UInt i = 0; i < size; ++i) {
-
-    if (surface_displacements(i) - shift - surface(i) < 0. &&    this->surface_tractions(i) < this->p_sat ){
-      shift = surface_displacements(i) - surface(i);
+        if (surface_displacements(i) - shift - surface(i) < 0. &&( this->p_sat < 0 || this->surface_tractions(i) < this->p_sat )){
+            shift = surface_displacements(i) - surface(i);
+        }
+    }
+#pragma omp parallel for
+    for (UInt i = 0; i < size; ++i) {
+        true_displacements(i) = surface_displacements(i) - shift;
     }
-  }
-
-  for (UInt i = 0; i < size; ++i) {
-  true_displacements(i) = surface_displacements(i) - shift;
- }
-  std::cout << std::scientific << std::setprecision(10)     << shift<<  std::fixed << std::endl;
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
diff --git a/src/maugis_adhesion_functional.cpp b/src/maugis_adhesion_functional.cpp
new file mode 100644
index 0000000..be2013a
--- /dev/null
+++ b/src/maugis_adhesion_functional.cpp
@@ -0,0 +1,104 @@
+/**
+ *
+
+ * @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 "maugis_adhesion_functional.hh"
+#include <cmath>
+
+/* -------------------------------------------------------------------------- */
+
+__BEGIN_TAMAAS__
+
+/* -------------------------------------------------------------------------- */
+
+MaugisAdhesionFunctional::MaugisAdhesionFunctional(BemFFTBase & bem):
+  Functional(bem),
+  term(bem.getSurface().size(), bem.getSurface().getL()),
+  gterm(bem.getSurface().size(), bem.getSurface().getL()),
+  surface_t(bem.getSurface().size(), bem.getSurface().getL())
+{}
+
+/* -------------------------------------------------------------------------- */
+
+MaugisAdhesionFunctional::~MaugisAdhesionFunctional()
+{}
+
+/* -------------------------------------------------------------------------- */
+
+Real MaugisAdhesionFunctional::computeF() {
+    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 (rho>=gap(i)){
+          this->term(i) = surface_energy *(1-gap(i)/rho);
+        }
+      }
+    }
+
+  return -SurfaceStatistics::computeSum(this->term);
+}
+
+/* -------------------------------------------------------------------------- */
+
+void MaugisAdhesionFunctional::computeGradFU() {
+    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 (rho>=gap(i)){
+          this->gterm(i) = surface_energy /rho;
+        }
+        else {this->gterm(i) =0;}
+      }
+    }
+
+  *this->gradF += this->gterm;
+}
+
+/* -------------------------------------------------------------------------- */
+
+void MaugisAdhesionFunctional::computeGradFP() {
+  SURFACE_FATAL("Adhesion functional (Maugis) cannot be differentiated with respect to P");
+
+}
+
+__END_TAMAAS__
diff --git a/python/bem.i b/src/maugis_adhesion_functional.hh
similarity index 50%
copy from python/bem.i
copy to src/maugis_adhesion_functional.hh
index 44fbea5..8fbf51e 100644
--- a/python/bem.i
+++ b/src/maugis_adhesion_functional.hh
@@ -1,50 +1,70 @@
 /**
  *
- * @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/>.
  *
  */
+
 /* -------------------------------------------------------------------------- */
 
+#ifndef __MAUGIS_ADHESION_FUNCTIONAL_HH__
+#define __MAUGIS_ADHESION_FUNCTIONAL_HH__
 
+#include "functional.hh"
+#include "bem_fft_base.hh"
 
-/* -------------------------------------------------------------------------- */
-/* Wrapping bem capabilities
-/* -------------------------------------------------------------------------- */
+__BEGIN_TAMAAS__
+
+/// Energy functional (in displacement) that contains the adhesion term from Maugis model
+class MaugisAdhesionFunctional : public Functional {
+public:
+  /// Constructor
+  explicit MaugisAdhesionFunctional(BemFFTBase & bem);
+
+  /// Destructor
+  virtual ~MaugisAdhesionFunctional();
+
+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();
+
+
+
+private:
+  ///  term in functional
+  Surface<Real> term;
+
+  ///  term in derivative functional
+  Surface<Real> gterm;
+
+  /// Surface tractions
+  Surface<Real> surface_t;
+};
+
+__END_TAMAAS__
 
-%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 "functional.hh"
-%include "meta_functional.hh"
-%include "elastic_energy_functional.hh"
-%include "complimentary_term_functional.hh"
-%include "exponential_adhesion_functional.hh"
-%include "fftransform.hh"
-
-namespace tamaas {
-%template(FFTransformReal) FFTransform<Real>;
-}
+#endif // __MAUGIS_ADHESION_FUNCTIONAL_HH__
diff --git a/tests/SConscript b/tests/SConscript
index bfeb4ef..2726165 100644
--- a/tests/SConscript
+++ b/tests/SConscript
@@ -1,20 +1,21 @@
 from os.path import join
 
 Import('main_env')
 
 test_files = Split("""
 run_tests.sh
 test_hertz_pressure.py
 test_westergaard.py
 test_fractal.py
 test_fractal.verified
 test_hertz_disp.py
+test_hertz_adhesion.py
 test_hertz_kato.py
 """)
 
 src_dir = "#/tests"
 build_dir = 'build-' + main_env['build_type'] + '/tests'
 
 for file in test_files:
     source = join(src_dir, file)
     Command(file, source, Copy("$TARGET", "$SOURCE"))
diff --git a/tests/test_hertz_adhesion.py b/tests/test_hertz_adhesion.py
new file mode 100644
index 0000000..6c5cb3b
--- /dev/null
+++ b/tests/test_hertz_adhesion.py
@@ -0,0 +1,134 @@
+# -*- coding: utf-8 -*-
+# -----------------------------------------------------------------------------
+# @author Lucas Frérot <lucas.frerot@epfl.ch>
+# @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/>.
+# -----------------------------------------------------------------------------
+
+import sys
+import tamaas as tm
+import numpy as np
+import scipy.optimize as scio 
+
+def constructHertzProfile(size, curvature):
+    radius = 1. / curvature
+    x = np.linspace(-0.5, 0.5, size)
+    y = np.linspace(-0.5, 0.5, size)
+    x, y = np.meshgrid(x, y)
+    surface=-(x**2 + y**2)/(2*radius)
+    return surface.copy()
+
+
+
+def main():
+    grid_size = 512
+    curvature = 0.5
+    effective_modulus = 1.
+    surface_energy = 0.0001
+    rho =  0.002
+	#Other quantities
+    radius=1./curvature
+    K=4*effective_modulus/3
+    KI=np.sqrt(surface_energy*2*effective_modulus)
+    tabor= surface_energy**(2/3)/(rho*effective_modulus**(2/3))*(1/curvature)**(1/3)
+    sig0=surface_energy/(rho)
+    lam=2*sig0/(np.pi*surface_energy*K**2/radius)**(1/3)
+
+
+    surface = constructHertzProfile(grid_size, curvature)
+    bem = tm.BemGigipol(surface)
+    bem.setDumpFreq(100)
+    functional = tm.MaugisAdhesionFunctional(bem)
+    functional.setParameter('rho', rho)
+    functional.setParameter('surface_energy', surface_energy)
+    bem.setEffectiveModulus(effective_modulus)
+    bem.addFunctional(functional)
+    bem.computeEquilibrium(1e-8, 0.00116480372339)
+
+    tractions = bem.getTractions()
+    true_gap = bem.getGap()
+    true_displacements = bem.getTrueDisplacements()
+
+    adhesive_pressure=np.zeros((grid_size,grid_size))
+    for i in np.arange(0,grid_size,1):
+        for j in np.arange(0,grid_size,1):
+            if (rho>true_gap[i,j]):
+                adhesive_pressure[i,j]=surface_energy/rho
+
+   # computed solution
+    p0=np.mean(tractions)
+    p_max=np.max(tractions)
+    contact_points=np.where(true_gap==0)
+    contact_area=np.size(contact_points)/2/(float(grid_size)**2)
+    ag=(contact_area/np.pi)**0.5
+    adh_points=np.where(adhesive_pressure!=0.)
+    adh_area=np.size(adh_points)/2/(float(grid_size)**2)
+    cg=(adh_area/np.pi)**0.5
+    delta=ag**2/radius-8*sig0/(3*K)*np.sqrt(cg**2-ag**2)
+	
+
+    # Exact solution
+    beta=(1-np.exp(lam/-0.924))/1.02
+    a_chap=1.54+((2.28*lam**1.3-1)/(2.28*lam**1.3+1))*0.279
+    L_chap=-7./4.+((4.04*lam**1.4-1)/(4.04*lam**1.4+1))/4.
+    a0=a_chap*(np.pi*surface_energy*radius**2/K)**(1./3.)
+    pc=L_chap*np.pi*surface_energy*radius
+    P=-0.00078
+    ac=a0*((beta+np.sqrt(1-P/pc))/(1+beta))**(2./3.)
+
+
+    def fonction1(m):
+        return np.sqrt(m**2-1)+m**2*np.arctan(np.sqrt(m**2-1))-np.pi*KI/(np.sqrt(np.pi*ac)*sig0)
+    m=scio.brenth(fonction1,1.,2.)
+    cc=m*ac
+
+
+    delta_cos=ac**2/radius-8*sig0/(3*K)*np.sqrt(cc**2-ac**2)
+
+
+    x = np.linspace(-0.5, 0.5, grid_size)
+    y = np.linspace(-0.5, 0.5, grid_size)
+    x, y = np.meshgrid(x, y)
+    a=ac
+    c=cc
+	# stress
+    sig_hertz=np.zeros((grid_size,grid_size))
+    disp=np.zeros((grid_size,grid_size))
+    for i in np.arange(0,grid_size,1):
+        for j in np.arange(0,grid_size,1):
+            if  a**2>(x[i,j]**2 + y[i,j]**2):
+                disp[i,j]=3*K/(2*np.pi*radius)*np.sqrt(a**2-(x[i,j]**2 + y[i,j]**2))-2*sig0/np.pi*np.arctan(np.sqrt((c**2-a**2)/(a**2-(x[i,j]**2 + y[i,j]**2))))
+            if  c**2>(x[i,j]**2 + y[i,j]**2)>a**2:
+                disp[i,j]=-sig0
+    sig_hertz=disp
+
+
+    a_error=abs((ac-ag)/ac)
+    c_error=abs((cc-cg)/cc)
+    delta_error=abs((delta-delta_cos)/delta)
+    pressure_error=abs((p0-np.mean(sig_hertz))/np.mean(sig_hertz))
+   
+
+    if a_error > 1.e-2 or c_error > 1.2e-2 or delta_error > 5.e-2 or pressure_error > 1.e-1 :
+        return 1
+    return 0
+if __name__ == "__main__":
+    sys.exit(main())