diff --git a/src/bem_gigi.cpp b/src/bem_gigi.cpp
index f5a483d..1994074 100644
--- a/src/bem_gigi.cpp
+++ b/src/bem_gigi.cpp
@@ -1,448 +1,441 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #include <vector>
 #include "surface.hh"
 #include "bem_gigi.hh"
 #include <iostream>
 #include <sstream>
 #include <fstream>
 #include <iomanip>
 #include <sstream>
 #include <cmath>
 /* -------------------------------------------------------------------------- */
 #define TIMER
 #include "surface_timer.hh"
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeEquilibrium(Real epsilon,
                      Real pressure,
                      int nthreads){
 
   this->setNumberOfThreads(nthreads);
   this->computeSpectralInfluenceOverDisplacement();
   this->surface_t = 0.;
   this->surface_r = 0.;
   this->surface_spectral_r = 0.;
   this->surface_w = 0.;
   this->surface_q = 0.;
   this->surface_spectral_q = 0.;
   this->pold = 0.;
   this->surface_displacements = 0.;
-  this->gradient_F = 0.;
   Real Gold = 1.;
   Real Rold = 1.;
   Real tau = 0.;
   Real f = 1e300;
   Real fPrevious = 1e300;
 
-  Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions);
+  // Real current_pressure = SurfaceStatistics::computeAverage(surface_tractions);
   surface_tractions = pressure;
   this->enforcePressureBalance(pressure);
 
   convergence_iterations.clear();
   nb_iterations = 0;
 
   while(f > epsilon && nb_iterations < max_iterations) {
      ++nb_iterations;
     fPrevious = f;
     this->computeDisplacements();
-    this->computeGaps();
-    Real gbar = this->computeMeanGapsInContact();
-    this->gap -= gbar;
-    //Juste pour le test
-    this->gradient_F = this->gap;
+    this->functional->computeGradF();
+    this->functional->centerGradFOnContactArea();
 
     Real R = this->computeR();
     this->updateW(R,Rold);
     Real alpha = this->computeAlpha();
     pold = this->surface_tractions;
     Rold = R;
     Gold = R;
     this->updatePressurea(alpha);
     tau=alpha;
     this->truncatePressure();
     UInt Iol=this->computeIol();
     //Inner loop
     while(Iol > 0 ) {
         ++nb_iterations;
         this->emptyoverlap(tau);
         this->enforcePressureBalance(pressure);
-        this->computeF();
-        f = this->getF();
+        this->functional->computeF();
+        f = this->functional->getF();
         fPrevious = f;
 
         this->computeDisplacements();
         this->computeGaps();
 
         Real gbar = this->computeMeanGapsInContact();
         this->gap -= gbar;
         Real G = this->computeG();
         this->updateT(G,Gold);
         Real tau = this->computeTau();
         pold = this->surface_tractions;
         Gold = G;
         this->updatePressureb(tau);
         this->truncatePressure();
         Iol= this->computeIol();
 
        }
 
     this->enforcePressureBalance(pressure);
-    this->computeF();
-    f = this->getF();
+    this->functional->computeF();
+    f = this->functional->getF();
 
     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::endl;
+      Real A =
+        SurfaceStatistics::computeContactAreaRatio(surface_tractions);
+      std::cout << std::scientific << std::setprecision(10)
+                << nb_iterations << " "
+                << f << " " << f-fPrevious << " " << A
+                <<  std::endl;
     }
 
     convergence_iterations.push_back(f);
 
 
   }
   return f;
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::truncatePressure(){
    STARTTIMER("truncatePressure");
 
    unsigned int n = surface.size();
    unsigned int size = n*n;
 #pragma omp parallel for
    for (unsigned int i = 0; i < size; ++i) {
 
     if (this->surface_tractions(i) < 0.) this->surface_tractions(i) = 0;
    }
      STOPTIMER("truncatePressure");
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 Real BemGigi::computeR(){
   STARTTIMER("computeR");
 
   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->surface_tractions(i) <= 0) continue;
-    Real val = this->gradient_F(i);
+    Real val = this->functional->getGradF()(i);
     res += val*val;
   }
 
   STOPTIMER("computeR");
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 Real BemGigi::computeIol(){
   STARTTIMER("computeIol");
 
   unsigned int n = surface.size();
   unsigned int size = n*n;
   UInt Iol = 0;
 #pragma omp parallel for
   for (unsigned int i = 0; i < size; ++i) {
         if (this->surface_tractions(i) <= 0. && this->gap(i) < 0.){
     ++Iol;}
   }
 
   STOPTIMER("Iol");
   return Iol;
 }
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::updateW(Real R, Real Rold){
   STARTTIMER("updateW");
 
   unsigned int n = surface.size();
   unsigned int size = n*n;
   Real factor = R/Rold;
 #pragma omp parallel for
   for (unsigned int i = 0; i < size; ++i) {
 
 
       if (this->surface_tractions(i) <= 0) this->surface_w(i) = 0.;
       else{
         this->surface_w(i) *= factor;
-        this->surface_w(i) += this->gradient_F(i);
+        this->surface_w(i) += this->functional->getGradF()(i);
       }
 
   }
 
   STOPTIMER("updateW");
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::updateT(Real G, Real Gold){
   STARTTIMER("updateT");
 
   unsigned int n = surface.size();
   unsigned int size = n*n;
   Real factor = G/Gold;
 #pragma omp parallel for
   for (unsigned int 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");
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeAlpha(){
   STARTTIMER("computeAlpha");
 
   surface_spectral_q = surface_w;
   surface_spectral_q.FFTTransform(nthreads);
 
   unsigned int n = surface.size();
   unsigned int size = n*n;
 #pragma omp parallel for
   for (unsigned int i = 0; i < size; ++i) {
     surface_spectral_q(i) *= this->surface_spectral_influence_disp(i);
   }
 
 
   surface_spectral_q.FFTITransform(surface_q,nthreads);
 
   Real qbar = 0;
   UInt nb_contact = 0;
 #pragma omp parallel for reduction(+: nb_contact,qbar)
   for (unsigned int i = 0; i < size; ++i) {
    if (this->surface_tractions(i) <= 0) continue;
     ++nb_contact;
     qbar += surface_q(i);
   }
   qbar /= nb_contact;
   surface_q -= qbar;
 
   Real alpha_sum1 = 0.,alpha_sum2 = 0.;
 
 #pragma omp parallel for reduction(+: alpha_sum1, alpha_sum2)
   for (unsigned int i = 0; i < size; ++i) {
-      if (this->surface_tractions(i) <= 0) continue;
-    alpha_sum1 += this->gradient_F(i)*surface_w(i);
+    if (this->surface_tractions(i) <= 0) continue;
+    alpha_sum1 += this->functional->getGradF()(i)*surface_w(i);
     alpha_sum2 += surface_q(i)*surface_w(i);
   }
   Real alpha = alpha_sum1/alpha_sum2;
 
   STOPTIMER("computeAlpha");
   return alpha;
 }
 
 
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::updatePressurea(Real alpha ){
    STARTTIMER("updatePressurea");
 
    unsigned int n = surface.size();
    unsigned int size = n*n;
 #pragma omp parallel for
    for (unsigned int i = 0; i < size; ++i) {
     this->surface_tractions(i) -= alpha*this->surface_w(i);
    }
 
    STOPTIMER("updatePressurea");
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::updatePressureb(Real tau ){
    STARTTIMER("updatePressureb");
 
    unsigned int n = surface.size();
    unsigned int size = n*n;
 #pragma omp parallel for
    for (unsigned int i = 0; i < size; ++i) {
     if (this->surface_tractions(i) <= 0) continue;
     this->surface_tractions(i) -= tau*this->surface_t(i);
    }
 
    STOPTIMER("updatePressureb");
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::emptyoverlap(Real tau){
    STARTTIMER("emptyoverlap");
 
    unsigned int n = surface.size();
    unsigned int size = n*n;
 #pragma omp parallel for
    for (unsigned int i = 0; i < size; ++i) {
-       if (this->surface_tractions(i) <= 0. && this->gap(i) < 0.){
-    this->surface_tractions(i) -= tau*this->gap(i);  }
+     if (this->surface_tractions(i) <= 0. && this->gap(i) < 0.){
+       this->surface_tractions(i) -= tau*this->gap(i);  }
    }
 
    STOPTIMER("emptyoverlap");
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::enforcePressureBalance(Real applied_pressure){
    STARTTIMER("enforcePressureBalance");
 
    unsigned int n = surface.size();
    unsigned int size = n*n;
    Real pressure = 0.;
 #pragma omp parallel for reduction(+: pressure)
    for (unsigned int i = 0; i < size; ++i) {
      pressure += this->surface_tractions(i);
    }
 
    pressure *= 1./size;
-   std::cerr << std::scientific << std::setprecision(10)
-            << "pressure " << pressure << std::endl;
+   // std::cerr << std::scientific << std::setprecision(10)
+            // << "pressure " << pressure << std::endl;
 #pragma omp parallel for
    for (unsigned int i = 0; i < size; ++i) {
      this->surface_tractions(i) *= applied_pressure/pressure;
    }
 
    STOPTIMER("enforcePressureBalance");
  }
 
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 
 void BemGigi::computeGaps() {
   UInt size = surface.size();
 
 #pragma omp parallel for
   for (UInt i=0; i < size*size; i++) {
     gap(i) = surface_displacements(i) - surface(i);
   }
 }
 
 
 /* -------------------------------------------------------------------------- */
 Real BemGigi::computeMeanGapsInContact(){
   STARTTIMER("computeMeanGapsInContact");
 
   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->surface_tractions(i) <= 0) continue;
     ++nb_contact;
     res += this->gap(i);
   }
 
   res /= nb_contact;
   STOPTIMER("computeMeanGapsInContact");
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeG(){
   STARTTIMER("computeG");
 
   unsigned int n = surface.size();
   unsigned int size = n*n;
   Real res = 0.;
 #pragma omp parallel for reduction(+: res)
   for (unsigned int i = 0; i < size; ++i) {
     if (this->surface_tractions(i) <= 0) continue;
     Real val = this->gap(i);
     res += val*val;
   }
 
   STOPTIMER("computeG");
   return res;
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 Real BemGigi::computeTau(){
   STARTTIMER("computeTau");
 
   surface_spectral_r = surface_t;
   surface_spectral_r.FFTTransform(nthreads);
 
   unsigned int n = surface.size();
   unsigned int size = n*n;
 #pragma omp parallel for
   for (unsigned int i = 0; i < size; ++i) {
     surface_spectral_r(i) *= this->surface_spectral_influence_disp(i);
   }
 
 
   surface_spectral_r.FFTITransform(surface_r,nthreads);
 
   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->surface_tractions(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->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;
 }
-
-
-
diff --git a/src/bem_gigi.hh b/src/bem_gigi.hh
index 93847af..fd6a02d 100644
--- a/src/bem_gigi.hh
+++ b/src/bem_gigi.hh
@@ -1,94 +1,91 @@
 #ifndef BEM_GIGI_H
 #define BEM_GIGI_H
 /* -------------------------------------------------------------------------- */
 #include "bem_fft_base.hh"
 /* -------------------------------------------------------------------------- */
 
 class BemGigi : public BemFFTBase {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
 
   BemGigi(Surface<Real> & p) :
     BemFFTBase(p),
     surface_t(p.size(),p.getL()),
     surface_r(p.size(),p.getL()),
     surface_spectral_r(p.size(),p.getL()),
     surface_w(p.size(),p.getL()),
     surface_q(p.size(),p.getL()),
-    gradient_F(p.size(),p.getL()),
     surface_spectral_q(p.size(),p.getL()),
     pold(p.size(),p.getL()){
     e_star = 1.;
     max_iterations = 10000;
   };
   virtual ~BemGigi(){};
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
   Real computeEquilibrium(Real epsilon,Real pressure,
               int nthreads = 1);
   // compute the non zero gap penalty functional
   Real computeR();
   Real computeG();
   // compute the conjugate step size
   Real computeAlpha();
   Real computeTau();
    // compute the overlapp space
   Real computeIol();
   void emptyoverlap(Real tau);
   // update the search direction
   void updateW(Real R, Real Rold);
    void updateT(Real G, Real Gold);
   // update the tractions
   void updatePressurea(Real alpha);
   void updatePressureb(Real tau);
 
   void truncatePressure();
   // enforce the applied and contact pressures balance
   void enforcePressureBalance(Real applied_pressure);
 
   Real computeMeanGapsInContact();
 
   void computeGaps();
 
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
 
   //! get the number of iterations
   UInt getNbIterations() const{return this->nb_iterations;};
   //! get the convergence steps of iterations
   const std::vector<Real> & getConvergenceIterations() const{return this->convergence_iterations;};
 
   Surface<Real> & getSurfaceT(){return surface_t;};
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
 
   //! exploration direction for the Gigi algo
   Surface<Real> surface_t;
   //! projected exploration direction for the Gigi algo
   Surface<Real> surface_r;
   //! projected exploration direction for the Gigi algo (fourier version)
   SurfaceComplex<Real> surface_spectral_r;
   //! exploration direction for the Gigi algo
   Surface<Real> surface_w;
   //! projected exploration direction for the Gigi algo
   Surface<Real> surface_q;
-  //! projected exploration direction for the Gigi algo
-  Surface<Real> gradient_F;
   //! projected exploration direction for the Gigi algo (fourier version)
   SurfaceComplex<Real> surface_spectral_q;
   //! previous pressure
   Surface<Real> pold;
 };
 
 #endif /* BEM_GIGI_H */
diff --git a/tests/run_tests.sh b/tests/run_tests.sh
index 4ac9a74..8c6f40b 100755
--- a/tests/run_tests.sh
+++ b/tests/run_tests.sh
@@ -1,25 +1,23 @@
 #!/bin/bash
 
-find .. -name "tamaas_environment.sh" -exec source '{}' +
-
 n_tests=0
 passed_tests=0
 for test_file in test_*.py; do
     ((n_tests++))
     echo "Testing $test_file"
     echo -n "------------ "
     python $test_file >/dev/null 2>&1
     result=$?
     if [ $result -ne 0 ]; then
         echo "failure"
     else
         echo "success"
         ((passed_tests++))
     fi
 done
 
 failed_tests=$((n_tests - passed_tests))
 echo ""
 echo "Test results :"
 echo "    - passed $passed_tests/$n_tests"
 echo "    - failed $failed_tests/$n_tests"
diff --git a/tests/test_hertzian.py b/tests/test_hertzian.py
index 6980a7d..9b29a0b 100644
--- a/tests/test_hertzian.py
+++ b/tests/test_hertzian.py
@@ -1,115 +1,115 @@
 #!/usr/bin/env python
 # coding: utf-8
 # -----------------------------------------------------------------------------
 # @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/>.
 # -----------------------------------------------------------------------------
 
 import sys
 import tamaas as tm
 import numpy as np
 import matplotlib.pyplot as plt
 
 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 = np.sqrt(radius**2 - x**2 - y**2)
     surface -= surface.mean()
     return surface.copy()
 
 def computeHertzDisplacement(e_star, contact_size, max_pressure, size):
     x = np.linspace(-0.5, 0.5, size)
     y = np.linspace(-0.5, 0.5, size)
     x, y = np.meshgrid(x, y)
     disp = np.pi * max_pressure / (4 * contact_size * e_star) * (2 * contact_size**2 - (x**2 + y**2))
     return disp.copy()
 
 def main():
     grid_size = 512
     curvature = 0.1
     effective_modulus = 1.
     load = 0.0001
 
     surface = constructHertzProfile(grid_size, curvature)
 
-    bem = tm.BemPolonski(surface)
+    bem = tm.BemGigi(surface)
     bem.setEffectiveModulus(effective_modulus)
     bem.computeEquilibrium(1e-12, load)
 
     tractions = bem.getTractions()
     displacements = bem.getDisplacements()
 
     # Testing contact area against Hertz solution for solids of revolution
     contact_area = tm.SurfaceStatistics.computeContactArea(tractions)
     hertz_contact_size = (3 * load / (4 * curvature * effective_modulus))**(1. / 3.)
     hertz_area = np.pi * hertz_contact_size**2
     area_error = np.abs(hertz_area - contact_area) / hertz_area
     print "Area error: {}".format(area_error)
 
     # Testing maximum pressure
     max_pressure = tractions.max()
     hertz_max_pressure = (6 * load * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi
     pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure
     print "Max pressure error: {}".format(pressure_error)
 
     # Testing displacements
     hertz_displacements = computeHertzDisplacement(effective_modulus,
                                                    hertz_contact_size,
                                                    hertz_max_pressure,
                                                    grid_size)
 
     # Selecing only the points that are in contact
     contact_indexes = [(i, j, tractions[i, j] > 0) for i in range(grid_size) for j in range(grid_size)]
     contact_indexes = map(lambda x: x[0:2], filter(lambda x: x[2], contact_indexes))
 
     # Displacements of bem are centered around the mean of the whole surface
     # and Hertz displacements are not centered, so we need to compute mean
     # on the contact zone for both arrays
     bem_mean = 0.
     hertz_mean = 0.
     for index in contact_indexes:
         bem_mean += displacements[index]
         hertz_mean += hertz_displacements[index]
 
     bem_mean /= len(contact_indexes)
     hertz_mean /= len(contact_indexes)
     # Correction applied when computing error
     correction = hertz_mean - bem_mean
 
     # Computation of error of displacement in contact zone
     error = 0.
     hertz_norm = 0.
     for index in contact_indexes:
         error += (hertz_displacements[index] - displacements[index] - correction)**2
         hertz_norm += (hertz_displacements[index] - hertz_mean)**2
 
     displacement_error = np.sqrt(error / hertz_norm)
     print "Displacement error (in contact zone): {}".format(displacement_error)
 
     if area_error > 1e-2 or pressure_error > 1e-2 or displacement_error > 1e-4:
         return 1
 
     return 0
 
 if __name__ == "__main__":
     sys.exit(main())