diff --git a/python/bem.i b/python/bem.i
index 2f58964..fe41438 100644
--- a/python/bem.i
+++ b/python/bem.i
@@ -1,53 +1,54 @@
 /**
  *
  * @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 "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"
 %include "fftransform_fftw.hh"
 
 namespace tamaas {
-%template(FFTransformFFTWReal) FFTransformFFTW<Real>;
+  %template(FFTransformReal) FFTransform<double,2>;
+  %template(FFTransformFFTWReal) FFTransformFFTW<double,2>;
 }
diff --git a/python/surface.i b/python/surface.i
index ccfce27..5407ca4 100644
--- a/python/surface.i
+++ b/python/surface.i
@@ -1,363 +1,366 @@
 /**
  *
  * @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 "bem_grid.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 "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>)
 
 /* -------------------------------------------------------------------------- */
 /* 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>;
   %template (Grid2dUInt) Grid<UInt, 2>;
 }
 
+%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 <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], Real L[d], UInt nb_components) {
+    void __init__(T * w, UInt n[d], UInt nb_components) {
       memcpy(this->n, n, d * sizeof(UInt));
-      memcpy(this->L, L, d * sizeof(Real));
+      memset(this->L, 1, d * sizeof(Real));
       this->nb_components = nb_components;
       this->data.wrapMemory(w, this->computeSize());
     }
 
-    GridForPython(PyObject * input, Real L[d] = {0}): Grid<T, d>() {
+    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
-      npy_intp grid_dims[d] = {0}; // tamaas dims
-      if (_n >= d)
-	memcpy(grid_dims, ndims, d * sizeof(npy_intp));
-      else {
+      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+1];
       else this->nb_components = 1;
 
       PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew(input);
       if (iter == NULL) throw;
-      this->__init__(( T *)(iter->dataptr), (UInt *)ndims, L, this->nb_components);
+      this->__init__(( T *)(iter->dataptr), grid_dims, this->nb_components);
       Py_DECREF(iter);
     }
 
-    ~GridForPython() {}
+    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);
   }
 }
 
 %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{
     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  )
 
 %grid_typemaps(tamaas::Grid2dReal, Real, NPY_DOUBLE, 2)
 %grid_typemaps(tamaas::Grid2dComplex, 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/bem_fft_base.cpp b/src/bem_fft_base.cpp
index 0bbee01..f92a785 100644
--- a/src/bem_fft_base.cpp
+++ b/src/bem_fft_base.cpp
@@ -1,270 +1,270 @@
 /**
  *
  * @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 "bem_fft_base.hh"
 /* -------------------------------------------------------------------------- */
 #define TIMER
 #include "surface_timer.hh"
 #include "surface_statistics.hh"
 #include "meta_functional.hh"
 #include "elastic_energy_functional.hh"
 #include "fft_plan_manager.hh"
 /* -------------------------------------------------------------------------- */
 #include <omp.h>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /* -------------------------------------------------------------------------- */
 
 BemFFTBase::BemFFTBase(Surface<Real> & p):
   BemInterface(p),
   functional(new MetaFunctional(*this)),
   surface_tractions(p.size(),p.getL()),
   surface_displacements(p.size(),p.getL()),
   surface_spectral_input(p.size(),p.getL()),
   surface_spectral_output(p.size(),p.getL()),
   surface_spectral_influence_disp(p.size(),p.getL()),
   true_displacements(p.size(),p.getL()),
   old_displacements(p.size(),p.getL()),
   gap(p.size(),p.getL()),
   e_star(std::nan("")),
   tmp_coeff(1./M_PI){
   //    this->setNumberOfThreads(this->nthreads);
   max_iterations = 5000;
   dump_freq = 100;
   this->functional->addFunctionalTerm(new ElasticEnergyFunctional(*this));
   std::cerr << this << " is setting up the surfaces";
   std::cerr << &p << " has size " << p.size() << std::endl;
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 BemFFTBase::~BemFFTBase() {
   // delete functional;
 }
 
 /* -------------------------------------------------------------------------- */
 
 const Surface<Real> & BemFFTBase::getTractions() const{
   return surface_tractions;
 }
 /* -------------------------------------------------------------------------- */
 const Surface<Real> & BemFFTBase::getDisplacements()const{
   return surface_displacements;
 }
 /* -------------------------------------------------------------------------- */
 const Surface<Real> & BemFFTBase::getSpectralInfluenceOverDisplacement(){
   return surface_spectral_influence_disp;
 }
 /* -------------------------------------------------------------------------- */
 void BemFFTBase::setEffectiveModulus(Real e_star){
   this->e_star = e_star;
   this->computeSpectralInfluenceOverDisplacement();
 }
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::computeSpectralInfluenceOverDisplacement(){
 
   STARTTIMER("computeSpectralInfluenceOverDisplacement");
 
   if (std::isnan(tmp_coeff)) SURFACE_FATAL("constant tmp_coeff is nan");
   if (std::isnan(e_star))    SURFACE_FATAL("constant e_star is nan");
 
   UInt n = surface.size();
 
   surface_spectral_influence_disp(0) = 0.;
 
   // everything but the external lines and central lines
 #pragma omp parallel for
   for (UInt i = 1; i < n/2; ++i) {
     for (UInt j = 1; j < n/2; ++j) {
 
       Real d = sqrt(i*i+j*j);
       surface_spectral_influence_disp(i,j)     = 1./d;
       surface_spectral_influence_disp(n-i,j)   = 1./d;
       surface_spectral_influence_disp(i,n-j)   = 1./d;
       surface_spectral_influence_disp(n-i,n-j) = 1./d;
 
     }
   }
   // external line (i or j == 0)
 #pragma omp parallel for
   for (UInt i = 1; i < n/2; ++i) {
     Real d = i;
     surface_spectral_influence_disp(i,0)   = 1./d;
     surface_spectral_influence_disp(n-i,0) = 1./d;
     surface_spectral_influence_disp(0,i)   = 1./d;
     surface_spectral_influence_disp(0,n-i) = 1./d;
 
   }
 
   //internal line (i or j == n/2)
 #pragma omp parallel for
   for (UInt i = 1; i < n/2; ++i) {
     Real d = sqrt(i*i+n/2*n/2);
     surface_spectral_influence_disp(i,n/2)   = 1./d;
     surface_spectral_influence_disp(n-i,n/2) = 1./d;
     surface_spectral_influence_disp(n/2,i)   = 1./d;
     surface_spectral_influence_disp(n/2,n-i) = 1./d;
   }
 
   {
     //i == n/2 j == n/2
     Real d = sqrt(2*n/2*n/2);
     surface_spectral_influence_disp(n/2,n/2) = 1./d;
     //i == 0 j == n/2
     d = n/2;
     surface_spectral_influence_disp(0,n/2) = 1./d;
     //i == n/2 j == 0
     surface_spectral_influence_disp(n/2,0) = 1./d;
   }
 
   UInt size = n*n;
   Real L = surface.getL();
 
  #pragma omp parallel for
    for (UInt i = 0; i < size; ++i) {
      surface_spectral_influence_disp(i) *= tmp_coeff/this->e_star*L;
    }
 
   STOPTIMER("computeSpectralInfluenceOverDisplacement");
 
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::computeDisplacementsFromTractions(){
 
   STARTTIMER("computeDisplacementsFromTractions");
 
   this->applyInfluenceFunctions(surface_tractions,surface_displacements);
 
   STOPTIMER("computeDisplacementFromTractions");
 
 }
 
 /* -------------------------------------------------------------------------- */
 void BemFFTBase::computeTractionsFromDisplacements(){
 
   STARTTIMER("computeDisplacementsFromTractions");
 
   this->applyInverseInfluenceFunctions(surface_displacements,surface_tractions);
 
   STOPTIMER("computeDisplacementFromTractions");
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::applyInfluenceFunctions(Surface<Real> & input, Surface<Real> & output){
 
   STARTTIMER("applyInfluenceFunctions");
 
-  FFTransform<Real> & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output);
+  FFTransform<Real, 2> & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output);
 
   plan1.forwardTransform();
 
   surface_spectral_output *= surface_spectral_influence_disp;
 
-  FFTransform<Real> & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output);
+  FFTransform<Real, 2> & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output);
 
   plan2.backwardTransform();
 
   STOPTIMER("applyInfluenceFunctions");
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::applyInverseInfluenceFunctions(Surface<Real> & input, Surface<Real> & output){
 
   STARTTIMER("applyInfluenceFunctions");
 
-  FFTransform<Real> & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output);
+  FFTransform<Real, 2> & plan1 = FFTPlanManager::get().createPlan(input,surface_spectral_output);
 
   plan1.forwardTransform();
 
   surface_spectral_output /= surface_spectral_influence_disp;
   surface_spectral_output(0) = 0;
 
-  FFTransform<Real> & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output);
+  FFTransform<Real, 2> & plan2 = FFTPlanManager::get().createPlan(output,surface_spectral_output);
 
   plan2.backwardTransform();
 
   STOPTIMER("applyInfluenceFunctions");
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::computeGaps() {
   UInt size = surface.size();
 
 #pragma omp parallel for
   for (UInt i = 0 ; i < size*size ; i++) {
     gap(i) = true_displacements(i) - surface(i);
   }
 }
 
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::computeTrueDisplacements(){
   this->true_displacements = this->surface_displacements;
 
   UInt n = surface.size();
   UInt size = n*n;
 
   Real shift = 1e300;
 
 
 
   for (UInt i = 0; i < size; ++i) {
     if (surface_displacements(i) - shift - surface(i) < 0.){
       shift = surface_displacements(i) - surface(i);
     }
   }
 
   for (UInt i = 0; i < size; ++i) {
     true_displacements(i) = surface_displacements(i) - shift;
   }
 
 
 
 }
 
 /* -------------------------------------------------------------------------- */
 
 void BemFFTBase::addFunctional(Functional *new_funct) {
   this->functional->addFunctionalTerm(new_funct);
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
diff --git a/src/bem_interface.hh b/src/bem_interface.hh
index 05ec6a4..215b89f 100644
--- a/src/bem_interface.hh
+++ b/src/bem_interface.hh
@@ -1,105 +1,105 @@
 /**
  *
  * @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/>.
  *
  */
 
 /* -------------------------------------------------------------------------- */
 #ifndef BEM_INTERFACE_H
 #define BEM_INTERFACE_H
 
 #include "surface.hh"
 #include "fftransform.hh"
 #include "fft_plan_manager.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 
 
 class BemInterface {
 
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
 
   BemInterface(Surface<Real> & surface):
     surface(surface),
     spectral_surface(surface.size(),surface.getL()) {
     this->computeSpectralSurface();
   };
   virtual ~BemInterface(){};
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 
 protected:
 
   //! compute the displacements from the tractions
   virtual void computeDisplacementsFromTractions(){SURFACE_FATAL("TO BE IMPLEMENTED")};
   //! compute the tractions from the displacements
   virtual void computeTractionsFromDisplacements(){SURFACE_FATAL("TO BE IMPLEMENTED")};
 
 public:
 
   //! compute the spectral surface
   virtual void computeSpectralSurface(){
-    FFTransform<Real> & plan = FFTPlanManager::get().createPlan(surface, spectral_surface);
+    FFTransform<Real, 2> & plan = FFTPlanManager::get().createPlan(surface, spectral_surface);
     plan.forwardTransform();
   };
 
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 
 public:
 
   //! return the actual tractions
   virtual const Surface<Real> & getTractions() const = 0;
   //! return the actual displacements
   virtual const Surface<Real> & getDisplacements() const = 0;
   //! return the manipulated surface
   const Surface<Real> & getSurface() const{return this->surface;};
   //! return the manipulated surface in fourier space
   const SurfaceComplex<Real> & getSpectralSurface() const{return this->spectral_surface;};
   //! return the manipulated surface
   void setSurface(Surface<Real> & surface){this->surface = surface;};
 
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 protected:
 
   //! the surface profile we work with
   Surface<Real> & surface;
   //! the surface profile in fourier space
   SurfaceComplex<Real> spectral_surface;
 
 };
 
 __END_TAMAAS__
 
 #endif /* BEM_INTERFACE_H */
diff --git a/src/fft_plan_manager.cpp b/src/fft_plan_manager.cpp
index c6adbb5..1a7c39d 100644
--- a/src/fft_plan_manager.cpp
+++ b/src/fft_plan_manager.cpp
@@ -1,140 +1,140 @@
 /**
  *
  * @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 "fft_plan_manager.hh"
 #include "fftransform_fftw.hh"
 #if defined(USE_CUFFT)
 #include "fftransform_cufft.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 
 
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 FFTPlanManager::FFTPlanManager(){
 }
 /* -------------------------------------------------------------------------- */
 
 FFTPlanManager::~FFTPlanManager(){
 }
 /* -------------------------------------------------------------------------- */
 
 
 FFTPlanManager & FFTPlanManager::get(){
   if (FFTPlanManager::singleton == NULL)
     FFTPlanManager::singleton = new FFTPlanManager();
 
   return *FFTPlanManager::singleton;
 }
 /* -------------------------------------------------------------------------- */
 
-FFTransform<Real> & FFTPlanManager::createPlan(Surface<Real> & input,
-                                               SurfaceComplex<Real> & output){
+FFTransform<Real, 2> & FFTPlanManager::createPlan(Grid<Real, 2> & input,
+						  GridHermitian<Real, 2> & output){
   auto index = std::make_pair(const_cast<Real*>(input.getInternalData()),
                               const_cast<complex*>(output.getInternalData()));
 
   auto it = plans.find(index);
   auto end = plans.end();
 
   if (it == end) {
 #if defined(USE_CUFFT)
-    plans[index] = new FFTransformCUFFT<Real>(input.size(),input,output);
+    plans[index] = new FFTransformCUFFT<Real, 2>(input.size(),input,output);
 #else
-    plans[index] = new FFTransformFFTW<Real>(input.size(),input,output);
+    plans[index] = new FFTransformFFTW<Real, 2>(input,output);
 #endif
   }
 
   return *plans[index];
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<>
 void FFTPlanManager::destroyPlan<Real>(const Real * some) {
   auto it = plans.begin();
   auto end = plans.end();
 
   for (; it != end ; ++it) {
     if (it->first.first == some) {
       delete it->second;
       plans.erase(it->first);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<>
 void FFTPlanManager::destroyPlan<complex>(const complex * some) {
   auto it = plans.begin();
   auto end = plans.end();
 
   for (; it != end ; ++it) {
     if (it->first.second == some) {
       delete it->second;
       plans.erase(it->first);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 /// FFT should not be used for any other type
 #define DESTROY(type) template<> void FFTPlanManager::destroyPlan<type>( \
                                  __attribute__((unused)) const type * some){}
 
 /// This list of types is defined in map_2d.cpp
 DESTROY(int);
 DESTROY(UInt);
 DESTROY(bool);
 DESTROY(unsigned long);
 
 #undef DESTROY
 
 /* -------------------------------------------------------------------------- */
 
-void FFTPlanManager::destroyPlan(Surface<Real> & input,
-                                 SurfaceComplex<Real> & output){
+void FFTPlanManager::destroyPlan(Grid<Real, 2> & input,
+                                 GridHermitian<Real, 2> & output){
 
   auto index = std::make_pair(const_cast<Real*>(input.getInternalData()),
 			      const_cast<complex*>(output.getInternalData()));
 
   auto it = plans.find(index);
   auto end = plans.end();
 
   if (it != end){
     delete (plans[index]);
     plans.erase(index);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 FFTPlanManager * FFTPlanManager::singleton = NULL;
 /* -------------------------------------------------------------------------- */
 
 
 
 __END_TAMAAS__
diff --git a/src/fft_plan_manager.hh b/src/fft_plan_manager.hh
index 053e4c5..ed8d35d 100644
--- a/src/fft_plan_manager.hh
+++ b/src/fft_plan_manager.hh
@@ -1,81 +1,81 @@
 /**
  *
  * @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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef FFT_PLAN_MANAGER_H
 #define FFT_PLAN_MANAGER_H
 /* -------------------------------------------------------------------------- */
 #include <map>
 #include "fftransform.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 class FFTPlanManager {
 
 
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 protected:
   FFTPlanManager();
 
 public:
   virtual ~FFTPlanManager();
 
 public:
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 
   /// Get singleton instance
   static FFTPlanManager & get();
 
   /// Create/retrieve a plan from two surfaces
-  FFTransform<Real> & createPlan(Surface<Real> & input, SurfaceComplex<Real> & output);
+  FFTransform<Real, 2> & createPlan(Grid<Real, 2> & input, GridHermitian<Real, 2> & output);
 
   /// Destroy a plan from two surfaces
-  void destroyPlan(Surface<Real> & input, SurfaceComplex<Real> & output);
+  void destroyPlan(Grid<Real, 2> & input, GridHermitian<Real, 2> & output);
 
   /// Destroy any plan containing a given data pointer
   template<typename T>
   void destroyPlan(const T * some);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
 
-  std::map<std::pair<Real*,complex*>,FFTransform<Real> *> plans;
+  std::map<std::pair<Real*,complex*>,FFTransform<Real, 2> *> plans;
 
 
   static FFTPlanManager * singleton;
 };
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 
 #endif /* FFT_PLAN_MANAGER_H */
diff --git a/src/fftransform.cpp b/src/fftransform.cpp
index b558e5e..51daf0b 100644
--- a/src/fftransform.cpp
+++ b/src/fftransform.cpp
@@ -1,60 +1,55 @@
 /**
  *
  * @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 "fftransform.hh"
 
 __BEGIN_TAMAAS__
 
-template<typename T>
-FFTransform<T>::FFTransform(UInt size,
-                            Surface<T> & real,
-                            SurfaceComplex<T> & spectral):
+template<typename T, UInt dim>
+FFTransform<T, dim>::FFTransform(Grid<T, dim> & real,
+				 GridHermitian<T, dim> & spectral):
   real(real),
-  spectral(spectral),
-  size(size),
-  reduced_size(size / 2 + 1)
+  spectral(spectral)
 {}
 
 /* -------------------------------------------------------------------------- */
 
-template<typename T>
-FFTransform<T>::~FFTransform() {}
+template<typename T, UInt dim>
+FFTransform<T, dim>::~FFTransform() {}
 
 /* -------------------------------------------------------------------------- */
 
-template<typename T>
-void FFTransform<T>::normalize() {
-#pragma omp parallel for
-  for (UInt i = 0 ; i < size * size ; ++i)
-    real(i) /= size * size;
+template<typename T, UInt dim>
+void FFTransform<T, dim>::normalize() {
+  real /= real.dataSize();
 }
 
 /* -------------------------------------------------------------------------- */
 
-template class FFTransform<Real>;
+template class FFTransform<Real, 2>;
 
 __END_TAMAAS__
diff --git a/src/fftransform.hh b/src/fftransform.hh
index 53ac91b..ac5164e 100644
--- a/src/fftransform.hh
+++ b/src/fftransform.hh
@@ -1,69 +1,66 @@
 /**
  *
  * @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 FFTRANSFORM_HH
 #define FFTRANSFORM_HH
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
-#include "surface.hh"
-#include "surface_complex.hh"
+#include "grid.hh"
+#include "grid_hermitian.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
-template <typename T>
+template <typename T, UInt dim>
 class FFTransform {
 public:
   /// Constructor
-  FFTransform(UInt size,
-              Surface<T> & real,
-              SurfaceComplex<T> & spectral);
+  FFTransform(Grid<T, dim> & real,
+              GridHermitian<T, dim> & spectral);
 
   /// Destructor
   virtual ~FFTransform();
 
 public:
   /// Perform FFT
   virtual void forwardTransform() = 0;
 
   /// Perform IFFT
   virtual void backwardTransform() = 0;
 
   /// Normalize the real surface after backward transform
   virtual void normalize();
 
 protected:
-  Surface<T> & real;
-  SurfaceComplex<T> & spectral;
-
-  UInt size, reduced_size;
+  Grid<T, dim> & real;
+  GridHermitian<T, dim> & spectral;
 };
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 #endif // FFTRANSFORM_HH
diff --git a/src/fftransform_fftw.cpp b/src/fftransform_fftw.cpp
index 3db360c..358a99d 100644
--- a/src/fftransform_fftw.cpp
+++ b/src/fftransform_fftw.cpp
@@ -1,79 +1,95 @@
 /**
  *
  * @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 "fftransform_fftw.hh"
 
 __BEGIN_TAMAAS__
 
-template<typename T>
-FFTransformFFTW<T>::FFTransformFFTW(UInt size,
-                                    Surface<T> & real,
-                                    SurfaceComplex<T> & spectral):
-  FFTransform<T>(size, real, spectral),
+template<typename T, UInt dim>
+FFTransformFFTW<T, dim>::FFTransformFFTW(Grid<T, dim> & real,
+					 GridHermitian<T, dim> & spectral):
+  FFTransform<T, dim>(real, spectral),
   forward_plan(NULL),
   backward_plan(NULL)
 {
+  TAMAAS_ASSERT(real.getNbComponents() == spectral.getNbComponents(),
+		"components for FFTW do not match");
+  // Data pointers
   Real * in = const_cast<Real *>(this->real.getInternalData());
   fftw_complex * out = reinterpret_cast<fftw_complex *>(this->spectral.getInternalData());
-  forward_plan = fftw_plan_dft_r2c_2d(this->size, this->size,
-                                      in,
-                                      out,
-                                      FFTW_ESTIMATE);
-  backward_plan = fftw_plan_dft_c2r_2d(this->size, this->size,
-                                       out,
-                                       in,
-                                       FFTW_ESTIMATE);
+
+  // fftw parameters
+  int rank = dim; // dimension of fft
+  const int * n = reinterpret_cast<const int *>(real.sizes()); // size of individual fft
+  int howmany = real.getNbComponents(); // how many fft to compute
+  int idist = 1, odist = 1; // components are next to each other in memory
+  int istride = real.getNbComponents(), ostride = real.getNbComponents();
+  int * inembeded = NULL, * onembeded = NULL; // row major
+
+
+  forward_plan = fftw_plan_many_dft_r2c(rank, n, howmany,
+					in, inembeded,
+					istride, idist,
+					out, onembeded,
+					ostride, odist,
+					FFTW_ESTIMATE);
+  backward_plan = fftw_plan_many_dft_c2r(rank, n, howmany,
+					 out, onembeded,
+					 ostride, odist,
+					 in, inembeded,
+					 istride, idist,
+					 FFTW_ESTIMATE);
 }
 
 /* -------------------------------------------------------------------------- */
 
-template<typename T>
-FFTransformFFTW<T>::~FFTransformFFTW() {
+template<typename T, UInt dim>
+FFTransformFFTW<T, dim>::~FFTransformFFTW() {
   fftw_destroy_plan(forward_plan);
   fftw_destroy_plan(backward_plan);
 }
 
 /* -------------------------------------------------------------------------- */
 
-template<typename T>
-void FFTransformFFTW<T>::forwardTransform() {
+template<typename T, UInt dim>
+void FFTransformFFTW<T, dim>::forwardTransform() {
   fftw_execute(forward_plan);
 }
 
 /* -------------------------------------------------------------------------- */
 
-template<typename T>
-void FFTransformFFTW<T>::backwardTransform() {
+template<typename T, UInt dim>
+void FFTransformFFTW<T, dim>::backwardTransform() {
   fftw_execute(backward_plan);
   this->normalize();
 }
 
 /* -------------------------------------------------------------------------- */
 
-template class FFTransformFFTW<Real>;
+template class FFTransformFFTW<Real, 2>;
 
 __END_TAMAAS__
diff --git a/src/fftransform_fftw.hh b/src/fftransform_fftw.hh
index 97cbe27..881cf82 100644
--- a/src/fftransform_fftw.hh
+++ b/src/fftransform_fftw.hh
@@ -1,61 +1,60 @@
 /**
  *
  * @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 FFTRANSFORM_FFTW_HH
 #define FFTRANSFORM_FFTW_HH
 
 #include "fftransform.hh"
 #include <fftw3.h>
 
 __BEGIN_TAMAAS__
 
-template <typename T>
-class FFTransformFFTW : public FFTransform<T> {
+template <typename T, UInt dim>
+class FFTransformFFTW : public FFTransform<T, dim> {
 public:
   /// Constructor
-  FFTransformFFTW(UInt size,
-                  Surface<T> & real,
-                  SurfaceComplex<T> & spectral);
+  FFTransformFFTW(Grid<T, dim> & real,
+		  GridHermitian<T, dim> & spectral);
 
   /// Destructor
   virtual ~FFTransformFFTW();
 
 public:
   /// Perform FFT
   void forwardTransform();
 
   /// Perform IFFT
   void backwardTransform();
 
 private:
   fftw_plan forward_plan;
   fftw_plan backward_plan;
 };
 
 __END_TAMAAS__
 
 #endif // FFTRANSFORM_FFTW_HH
diff --git a/src/grid_hermitian.hh b/src/grid_hermitian.hh
index d99a632..5982cfe 100644
--- a/src/grid_hermitian.hh
+++ b/src/grid_hermitian.hh
@@ -1,163 +1,161 @@
 /**
  *
  * @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 __GRID_HERMITIAN_HH__
 #define __GRID_HERMITIAN_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid.hh"
 #include <complex>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 static std::complex<Real> dummy;
 
 /**
  * @brief Multi-dimensional, multi-component herimitian array
  */
 template <typename T, UInt dim>
 class GridHermitian : public Grid<std::complex<T>, dim> {
 /* -------------------------------------------------------------------------- */
 /* Constructors                                                               */
 /* -------------------------------------------------------------------------- */
 public:
   /// Constructor by default (empty array)
   GridHermitian();
   /// Constructor
   GridHermitian(const UInt n[dim], const Real L[dim], UInt nb_components);
   /// Destructor
   virtual ~GridHermitian();
 
 /* -------------------------------------------------------------------------- */
 /* Access operators                                                           */
 /* -------------------------------------------------------------------------- */
   template <typename... T1>
   inline std::complex<T> & operator()(T1... args);
   template <typename... T1>
   inline const std::complex<T> & operator()(T1... args) const;
 
 private:
   template <typename... T1>
   inline void packTuple(UInt * t, UInt i, T1... args) const;
   template <typename... T1>
   inline void packTuple(UInt * t, UInt i) const;
 
 };
 
 /* -------------------------------------------------------------------------- */
 /* Inline function definitions                                                */
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline void GridHermitian<T, dim>::packTuple(UInt * t, UInt i, T1... args) const {
   *t = i;
   packTuple(t+1, args...);
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline void GridHermitian<T, dim>::packTuple(UInt * t, UInt i) const {
   *t = i;
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline std::complex<T> & GridHermitian<T, dim>::operator()(T1... args) {
   TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1,
 		"operator() does not match dimension");
   if (sizeof...(T1) == 1 && dim != 1) {
     UInt tuple[1] = {0};
     packTuple(tuple, args...);
     if (tuple[0] >= this->data.size()) {
       TAMAAS_DEBUG_EXCEPTION("out of bonds access on hermitian surface");
       dummy = std::complex<T>(0,0);
       return dummy;
     }
     else return this->Grid<std::complex<T>, dim>::operator()(args...);
   }
 
   UInt tuple[dim+1] = {0};
   packTuple(tuple, args...);
 
   // reverse tuple if necessary
   if (tuple[dim-1] >= this->n[dim-1]) {
     for (UInt i = 0 ; i < dim ; i++) {
       if (tuple[i] && i != dim-1)
 	tuple[i] = this->n[i]-tuple[i];
       else if (tuple[i] && i == dim-1)
 	tuple[i] = 2*(this->n[i]-1) - tuple[i];
     }
 
     dummy = std::conj(this->Grid<std::complex<T>, dim>::operator()(tuple));
     return dummy;
   } else {
     return this->Grid<std::complex<T>, dim>::operator()(tuple);
   }
 }
 
 template <typename T, UInt dim>
 template <typename... T1>
 inline const std::complex<T> & GridHermitian<T, dim>::operator()(T1... args) const {
   TAMAAS_ASSERT(sizeof...(T1) == dim + 1 || sizeof...(T1) == 1,
 		"operator() does not match dimension");
   if (sizeof...(T1) == 1 && dim != 1) {
     UInt tuple[1] = {0};
     packTuple(tuple, args...);
     if (tuple[0] >= this->data.size()) {
       TAMAAS_DEBUG_EXCEPTION("out of bonds access on hermitian surface");
       dummy = std::complex<T>(0,0);
       return dummy;
     }
     else return this->Grid<std::complex<T>, dim>::operator()(args...);
   }
 
-  UInt tuple[dim+1] = {0}, reverse_tuple[dim+1] = {0};
+  UInt tuple[dim+1] = {0};
   packTuple(tuple, args...);
 
   // reverse tuple if necessary
   if (tuple[dim-1] >= this->n[dim-1]) {
     for (UInt i = 0 ; i < dim ; i++) {
       if (tuple[i] && i != dim-1)
-	reverse_tuple[i] = this->n[i]-tuple[i];
+	tuple[i] = this->n[i]-tuple[i];
       else if (tuple[i] && i == dim-1)
-	reverse_tuple[i] = 2*(this->n[i]-1) - tuple[i];
+	tuple[i] = 2*(this->n[i]-1) - tuple[i];
     }
+
+    dummy = std::conj(this->Grid<std::complex<T>, dim>::operator()(tuple));
+    return dummy;
   } else {
-    for (UInt i = 0 ; i < dim ; i++)
-      reverse_tuple[i] = tuple[i];
+    return this->Grid<std::complex<T>, dim>::operator()(tuple);
   }
-
-  dummy = std::conj(this->Grid<std::complex<T>, dim>::operator()(reverse_tuple));
-  return dummy;
 }
 
-
 __END_TAMAAS__
 
 
 
 #endif // __GRID_HERMITIAN_HH__
diff --git a/src/surface_generator_filter.cpp b/src/surface_generator_filter.cpp
index 8dbe39a..4dd4514 100644
--- a/src/surface_generator_filter.cpp
+++ b/src/surface_generator_filter.cpp
@@ -1,147 +1,147 @@
 /**
  *
  * @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> & plan = FFTPlanManager::get().createPlan(*surface, *surface_spectral);
+  FFTransform<Real, 2> & plan = FFTPlanManager::get().createPlan(*surface, *surface_spectral);
   plan.backwardTransform();
 }
 /* -------------------------------------------------------------------------- */
 
 
 void SurfaceGeneratorFilter::generateWhiteNoiseFFT(){
   white_noise_FFT = new SurfaceComplex<Real>(sources->size(), sources->getL());
-  FFTransform<Real> & plan = FFTPlanManager::get().createPlan(*sources, *white_noise_FFT);
+  FFTransform<Real, 2> & plan = FFTPlanManager::get().createPlan(*sources, *white_noise_FFT);
   plan.forwardTransform();
   // white_noise_FFT->dumpToParaview("white noise FFT",NO_REPLICATION);
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 void SurfaceGeneratorFilter::generateHFFT(){
-  FFTransform<Real> & plan = FFTPlanManager::get().createPlan(*h_coeff, *h_coeff_power);
+  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();
 }
 /* -------------------------------------------------------------------------- */
 
 
 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__
diff --git a/src/surface_statistics.hh b/src/surface_statistics.hh
index 39c22b4..ca53377 100644
--- a/src/surface_statistics.hh
+++ b/src/surface_statistics.hh
@@ -1,721 +1,721 @@
 /**
  *
  * @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/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 #ifndef __SURFACE_STATISTICS_HH__
 #define __SURFACE_STATISTICS_HH__
 /* -------------------------------------------------------------------------- */
 #include "surface.hh"
 #include "fftransform.hh"
 #include "fft_plan_manager.hh"
 #include <fftw3.h>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 
 class SurfaceStatistics {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
 
   SurfaceStatistics(){};
   virtual ~SurfaceStatistics(){};
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
 
   inline static Real computeMaximum(const Surface<Real> & s);
   inline static Real computeMinimum(const Surface<Real> & s);
   inline static UInt computeAverage(const Surface<UInt> & s,
 				    UInt epsilon=0){SURFACE_FATAL("TO IMPLEMENT");};
 
   inline static unsigned long computeAverage(const Surface<unsigned long> & s,
 					     unsigned long epsilon=0){SURFACE_FATAL("TO IMPLEMENT");};
 
 
   inline static Real computeAverage(const Surface<Real> & s,
 				    Real epsilon= -1.*std::numeric_limits<Real>::max()){
     return computeAverage<Real>(s,epsilon);
   }
 
   template <typename T>
   inline static T computeAverage(const Surface<T> & s,
 				 T epsilon= -1.*std::numeric_limits<T>::max());
   inline static Real computeStdev(const Surface<Real> & s, Real epsilon=-1. * std::numeric_limits<Real>::max());
 
   inline static Real computeSkewness(const Surface<Real> & surface, Real avg, Real std);
   inline static Real computeKurtosis(const Surface<Real> & surface, Real avg, Real std);
   inline static Real computeRMSSlope(const Surface<Real> & surface);
   inline static Real computeSpectralRMSSlope(Surface<Real> & surface);
   inline static Real computeSpectralMeanCurvature(Surface<Real> & surface);
   inline static Real computeSpectralStdev(Surface<Real> & surface);
   //  inline static void computeStatistics(Surface<Real> & surface);
   inline static std::vector<Real> computeMoments(Surface<Real> & surface);
   inline static Real computeAnalyticFractalMoment(UInt order,UInt k1, UInt k2, Real hurst, Real rms, Real L);
   //! compute perimeter
   inline static Real computePerimeter(const Surface<Real> & surface, Real eps = 1e-10);
   //! get real contact area
   inline static Real computeContactArea(const Surface<Real> & surface);
   //! get real contact area
   inline static Real computeContactAreaRatio(const Surface<Real> & surface);
 
 
   template <bool consider_zeros=false>
   inline static std::vector<int>
   computeDistribution(const Surface<Real> & surface,
 		      unsigned int n_bins,
 		      Real min_value=std::numeric_limits<Real>::max(),
 		      Real max_value=-1*std::numeric_limits<Real>::max());
 
   inline static std::vector<complex>
   computeSpectralDistribution(const Surface<Real> & surface,
 			      unsigned int n_bins,
 			      UInt cutoff=std::numeric_limits<UInt>::max());
 
 
 
   //! compute average amplitude
   inline static Real computeAverageAmplitude(const Surface<Real> & surface);
   inline static Real computeAverageAmplitude(const SurfaceComplex<Real> & surface);
 
 
   //! sum all heights
   inline static Real computeSum(const Surface<Real> & surface);
   //! compute auto correlation
   inline static Surface<Real> computeAutocorrelation(Surface<Real> & surface);
   //! compute powerspectrum
   inline static SurfaceComplex<Real> computePowerSpectrum(Surface<Real> & surface);
 
 
 };
 
 
 /* -------------------------------------------------------------------------- */
 /* inline functions                                                           */
 /* -------------------------------------------------------------------------- */
 
 //#include "surface_statistics_inline_impl.cc"
 
 /* -------------------------------------------------------------------------- */
 
 template <bool consider_zeros>
 std::vector<int>
 SurfaceStatistics::computeDistribution(const Surface<Real> & surface,
 				       unsigned int n_bins,
 				       Real min_value,Real max_value){
 
   if (n_bins == 0) SURFACE_FATAL("badly set number of bins");
 
   Real min_heights;
   Real max_heights;
 
   std::vector<int> bins;
   bins.resize(n_bins);
   bins.assign(n_bins,0);
 
 
   min_heights = min_value;
   max_heights = max_value;
 
   if (max_value == -1.*std::numeric_limits<Real>::max() )
     max_heights = computeMaximum(surface);
   if (min_value == std::numeric_limits<Real>::max() )
     min_heights = computeMinimum(surface);
 
   int counter = 0;
   const int n = surface.size();
   for (int i = 0 ; i < n*n ; ++i)
     {
       Real value = surface(i);
       if (!consider_zeros && value < zero_threshold) continue;
       if (value > max_heights || value < min_heights)
 	SURFACE_FATAL("out of bound value: " << min_heights << " <? " << value
 		      << " <? " << max_heights << " : bounds probably badly chosen");
       Real ind = (value - min_heights) / (max_heights - min_heights);
       unsigned int index = ind*n_bins;
       if (index == n_bins) --index;
       ++bins[index];
       ++counter;
     }
   if (consider_zeros && counter != n*n)
     SURFACE_FATAL("problem in counting points");
 
   //  Real offset = (max_heights - min_heights)/n_bins;
   //Real x = 0;
   //  Real f1 = 0, f2 = 0;
   //  for (unsigned int i = 0 ; i < n_bins ; ++i)
   //  {
   //    x  = i*offset+min_heights;
   //f1 = bins[i]/Real(n*n);
   //f2 = bins[i]/Real(counter);
   //}
   Real min_non_zero = std::numeric_limits<Real>::max();
 #pragma omp parallel for reduction(min: min_non_zero)
   for (int i = 0 ; i < n*n ; ++i) {
     Real value = surface(i);
     if (value < min_non_zero && value > 1e-14)
       min_non_zero = value;
   }
   return bins;
 }
 /* -------------------------------------------------------------------------- */
 std::vector<complex>
 SurfaceStatistics::computeSpectralDistribution(const Surface<Real> & surface,
 					       unsigned int n_bins,
 					       UInt cutoff){
 
   if (n_bins == 0) SURFACE_FATAL("badly set number of bins");
 
 
   const int n = surface.size();
 
 
   Real min_value = SurfaceStatistics::computeMinimum(surface);
   Real max_value = SurfaceStatistics::computeMaximum(surface);
 
   std::vector<complex> distrib_fft(n_bins);
   const int N = n_bins;
   Real L = max_value - min_value;
 
 #pragma omp parallel for
   for (int k = 0 ; k < N ; ++k) distrib_fft[k] = 0.;
 
 #pragma omp parallel for
   for (int k = 1 ; k < N/2 ; ++k) {
     for (int i = 0 ; i < n*n ; ++i) {
       Real value = surface(i);
       Real angle = -2.*M_PI*k*(value-min_value)/L;
       complex phase = std::polar(1.,angle);
       distrib_fft[k]   += phase;
       distrib_fft[N-k] += phase;
     }
   }
 
   for (int i = 0 ; i < n*n ; ++i) {
     Real value = surface(i);
     Real angle = -2.*M_PI*N/2*(value-min_value)/L;
     complex phase = std::polar(1.,angle);
     distrib_fft[N/2]   += phase;
   }
 
 #pragma omp parallel for
   for (int k = 1 ; k < N ; ++k) distrib_fft[k] /= n*n;
 
   distrib_fft[0] = std::polar(1.,0.);
 
 #pragma omp parallel for
   for (int k = 1 ; k < N/2 ; ++k) {
     if (k >= (int)cutoff) {
       distrib_fft[k] = 0.;
       distrib_fft[N-k] = 0.;
     }
   }
 
   //proceed to inverse fourier transform
 
 #pragma omp parallel for
   for (int k = 0 ; k < N ; ++k) {
     distrib_fft[k] = std::conj(distrib_fft[k]);
   }
 
   fftw_complex * in  = (fftw_complex*)&distrib_fft[0];
   fftw_plan p = fftw_plan_dft_1d(N, in, in ,FFTW_FORWARD, FFTW_ESTIMATE);
   fftw_execute(p);
   fftw_destroy_plan(p);
 
   std::vector<complex> distrib;
   distrib.resize(N);
 
 #pragma omp parallel for
   for (int k = 0 ; k < N ; ++k) {
     distrib[k] = std::conj(distrib_fft[k]);
   }
 
   return distrib;
 
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeAverageAmplitude(const Surface<Real> & surface){
   Real average_amplitude = 0;
   UInt n = surface.size();
 #pragma omp parallel for collapse(2), reduction(+: average_amplitude)
   for (UInt i = 0 ; i < n ; ++i)
     for (UInt j = 0 ; j < n ; ++j)
       average_amplitude += std::abs(surface(i,j));
 
   average_amplitude *= 1./n/n;
   return average_amplitude;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeAverageAmplitude(const SurfaceComplex<Real> & surface){
   Real average_amplitude = 0;
   UInt n = surface.size();
 #pragma omp parallel for collapse(2), reduction(+: average_amplitude)
   for (UInt i = 0 ; i < n ; ++i)
     for (UInt j = 0 ; j < n ; ++j)
       average_amplitude += std::abs(surface(i,j));
 
   average_amplitude *= 1./n/n;
   return average_amplitude;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeContactArea(const Surface<Real> & surface) {
   Real L = surface.getL();
   return SurfaceStatistics::computeContactAreaRatio(surface)*L*L;
 }
 
 
 /* -------------------------------------------------------------------------- */
 Real SurfaceStatistics::computeContactAreaRatio(const Surface<Real> & surface) {
   UInt area = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for collapse(2), reduction(+:area)
   for (UInt i = 0; i<n; i++) {
     for (UInt j = 0; j<n; j++) {
       if (surface (i,j) > 1e-10)
         area++;
     }
   }
   return area/Real(n*n);
 }
 
 
 
 
 /* -------------------------------------------------------------------------- */
 
 Surface<Real> SurfaceStatistics::computeAutocorrelation(Surface<Real> & surface){
   SurfaceComplex<Real> correlation_heights = computePowerSpectrum(surface);
   Surface<Real> acf(correlation_heights.size(), correlation_heights.getL());
-  FFTransform<Real> & plan = FFTPlanManager::get().createPlan(acf, correlation_heights);
+  FFTransform<Real, 2> & plan = FFTPlanManager::get().createPlan(acf, correlation_heights);
   plan.backwardTransform();
   return acf;
 }
 /* -------------------------------------------------------------------------- */
 
 
 SurfaceComplex<Real> SurfaceStatistics::computePowerSpectrum(Surface<Real> & surface){
   SurfaceComplex<Real> power_spectrum(surface.size(), surface.getL());
-  FFTransform<Real> & plan = FFTPlanManager::get().createPlan(surface, power_spectrum);
+  FFTransform<Real, 2> & plan = FFTPlanManager::get().createPlan(surface, power_spectrum);
   plan.forwardTransform();
 
   UInt n = power_spectrum.size();
   Real factor = 1./(n*n);
 
 #pragma omp parallel for collapse(2)
   for (UInt i = 0 ; i < n ; ++i)
     for (UInt j = 0 ; j < n ; ++j)
       power_spectrum(i,j) *= factor;
 
   power_spectrum.makeItRealBySquare();
   return power_spectrum;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 T SurfaceStatistics::computeAverage(const Surface<T> & surface, T epsilon){
 
   Real avg = 0.0;
   UInt cpt = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+: avg, cpt)
   for (UInt i = 0; i < n*n; i++) {
     Real val = surface(i);
     if (std::abs(val) > epsilon){
       avg += val;
       ++cpt;
     }
   }
   avg /= cpt;
   return avg;
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeSkewness(const Surface<Real> & surface,
                                         Real avg, Real std){
   Real skw = 0.0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+: skw)
   for (UInt i = 0; i < n*n; i++) {
     Real tmp = surface(i) - avg;
     skw += (tmp*tmp*tmp);
   }
   skw /= (n*std*std*std);
   return skw;
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeKurtosis(const Surface<Real> & surface,
                                         Real avg, Real std){
 
   Real kur=0.0;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+: kur)
   for (UInt i = 0; i < n*n; i++) {
     Real tmp = surface.size() - avg;
     kur += (tmp*tmp*tmp*tmp);
   }
   kur /= (n*std*std*std*std);
   return kur-3;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeMaximum(const Surface<Real> & surface){
   Real _max = -1. * std::numeric_limits<Real>::max();
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(max: _max)
   for (UInt i = 0; i < n*n; i++) {
     _max = std::max(surface(i),_max);
   }
   return _max;
 }
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeMinimum(const Surface<Real> & surface){
   Real _min = std::numeric_limits<Real>::max();
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(min: _min)
   for (UInt i = 0; i < n*n; i++) {
     _min = std::min(surface(i),_min);
   }
   return _min;
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeRMSSlope(const Surface<Real> & surface){
 
   UInt n = surface.size();
   Surface<Real> grad(n,1.);
 
   for (UInt i = 0, m = 0; i < n; i++) {
     for (UInt j = 0; j < n; j++) {
       UInt k = i+n*j;
 
       //      xl = ((i==0)?(k+n-1):(k-1));
       UInt xh = ((i==(n-1))?(k-n+1):(k+1));
       //      yl = ((j==0)?(k+n*n-n):(k-n));
       UInt yh = ((j==(n-1))?(k-n*n+n):(k+n));
 
 
       Real gx = surface(xh) - surface(k);
       Real gy = surface(yh) - surface(k);
       grad(m++) = gx*gx + gy*gy;
     }
   }
 
   Real avg = SurfaceStatistics::computeAverage(grad);
   Real slp = sqrt(avg); /* According to paper Hyun PRE 70:026117 (2004) */
   Real scale_factor = surface.getL()/n;
   return slp/scale_factor;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSpectralRMSSlope(Surface<Real> & surface){
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   UInt n = surface.size();
 
   Real rms_slopes = 0;
 #pragma omp parallel for collapse(2), reduction(+:rms_slopes)
   for (UInt i = 1 ; i < n/2 ; ++i)
     for (UInt j = 1 ; j < n/2 ; ++j){
       rms_slopes += 1.*((i*i)+(j*j))*power(i,j).real();
       rms_slopes += 1.*((i*i)+(j*j))*power(n-i,n-j).real();
       rms_slopes += 1.*((i*i)+(j*j))*power(n-i,j).real();
       rms_slopes += 1.*((i*i)+(j*j))*power(i,n-j).real();
     }
 
   // external line (i or j == 0)
 #pragma omp parallel for reduction(+:rms_slopes)
   for (UInt i = 1 ; i < n/2 ; ++i){
     rms_slopes += 1.*(i*i)*power(i,0).real();
     rms_slopes += 1.*(i*i)*power(n-i,0).real();
     rms_slopes += 1.*(i*i)*power(0,i).real();
     rms_slopes += 1.*(i*i)*power(0,n-i).real();
   }
 
   //internal line (i or j == n/2)
 #pragma omp parallel for reduction(+:rms_slopes)
   for (UInt i = 1; i < n/2; ++i) {
     rms_slopes += 0.5*(i*i+n*n/4)*power(i,n/2).real();
     rms_slopes += 0.5*(i*i+n*n/4)*power(n-i,n/2).real();
     rms_slopes += 0.5*(i*i+n*n/4)*power(n/2,i).real();
     rms_slopes += 0.5*(i*i+n*n/4)*power(n/2,n-i).real();
   }
 
   {
     //i == n/2 j == n/2
     rms_slopes += 0.25*(n*n/2)*power(n/2,n/2).real();
     //i == 0 j == n/2
     rms_slopes += 1.*(n*n/4)*power(0,n/2).real();
     //i == n/2 j == 0
     rms_slopes += 1.*(n*n/4)*power(n/2,0).real();
   }
   rms_slopes = 2.*M_PI/n*sqrt(rms_slopes);
   Real scale_factor = surface.getL()/n;
   return rms_slopes/scale_factor;
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeStdev(const Surface<Real> & surface, Real epsilon){
   UInt n = surface.size();
   UInt cpt = 0;
   Real std = 0.0;
   Real avg = 0.0;
 
 #pragma omp parallel for reduction(+:avg, std, cpt)
   for (UInt i = 0; i < n*n; i++) {
     Real val = surface(i);
     if (val > epsilon){
       avg += val;
       std += val*val;
       ++cpt;
     }
   }
   avg /= cpt;
   return sqrt(std/cpt - avg*avg);
 }
 /* -------------------------------------------------------------------------- */
 
 
 Real SurfaceStatistics::computeSpectralStdev(Surface<Real> & surface){
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   Real rms_heights = 0;
   UInt n = surface.size();
 
 #pragma omp parallel for collapse(2), reduction(+:rms_heights)
   for (UInt i = 0 ; i < n ; ++i)
     for (UInt j = 0 ; j < n ; ++j){
       rms_heights += power(i,j).real();
     }
   rms_heights = sqrt(rms_heights);
   return rms_heights;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeSpectralMeanCurvature(Surface<Real> & surface){
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   UInt n = surface.size();
 
   Real rms_curvatures = 0;
 #pragma omp parallel for collapse(2), reduction(+:rms_curvatures)
   for (UInt i = 1 ; i < n/2 ; ++i)
     for (UInt j = 1 ; j < n/2 ; ++j){
       Real coeff = ((i*i)+(j*j));
       coeff *= coeff;
       rms_curvatures += 1.*coeff*power(i,j).real();
       rms_curvatures += 1.*coeff*power(n-i,n-j).real();
       rms_curvatures += 1.*coeff*power(n-i,j).real();
       rms_curvatures += 1.*coeff*power(i,n-j).real();
     }
 
   // external line (i or j == 0)
 #pragma omp parallel for reduction(+:rms_curvatures)
   for (UInt i = 1 ; i < n/2 ; ++i){
     Real coeff = i*i;
     coeff *= coeff;
     rms_curvatures += 4./3.*coeff*power(i,0).real();
     rms_curvatures += 4./3.*coeff*power(n-i,0).real();
     rms_curvatures += 4./3.*coeff*power(0,i).real();
     rms_curvatures += 4./3.*coeff*power(0,n-i).real();
   }
 
   //internal line (i or j == n/2)
 #pragma omp parallel for reduction(+:rms_curvatures)
   for (UInt i = 1; i < n/2; ++i) {
     Real coeff = (i*i+n*n/4);
     coeff *= coeff;
     rms_curvatures += 0.5*coeff*power(i,n/2).real();
     rms_curvatures += 0.5*coeff*power(n-i,n/2).real();
     rms_curvatures += 0.5*coeff*power(n/2,i).real();
     rms_curvatures += 0.5*coeff*power(n/2,n-i).real();
   }
 
   {
     //i == n/2 j == n/2
     rms_curvatures += 0.25*(n*n/2)*(n*n/2)*power(n/2,n/2).real();
     //i == 0 j == n/2
     rms_curvatures += 2./3.*(n*n/4)*(n*n/4)*power(0,n/2).real();
     //i == n/2 j == 0
     rms_curvatures += 2./3.*(n*n/4)*(n*n/4)*power(n/2,0).real();
   }
 
   rms_curvatures *= 9./4.;
   rms_curvatures = M_PI/n/n * sqrt(rms_curvatures);
 
   return rms_curvatures;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computePerimeter(const Surface<Real> & surface, Real eps) {
   int perimeter = 0;
   UInt n = surface.size();
 
 
 #pragma omp parallel for collapse(2), reduction(+:perimeter)
   for (UInt i = 0 ; i < n ; ++i) {
     for (UInt j = 0 ; j < n ; ++j) {
       if ( j > 0 && ( ( surface(i,(j-1)) < eps && surface(i,j) > eps ) ||
                       ( surface(i,(j-1)) > eps && surface(i,j) < eps ) ) )
         perimeter++;
       if ( i > 0 && ( ( surface(j,(i-1)) < eps && surface(j,i) > eps ) ||
                       ( surface(j,(i-1)) > eps && surface(j,i) < eps )) )
         perimeter++;
     }
   }
 
   return perimeter/Real(n);
 }
 /* -------------------------------------------------------------------------- */
 Real SurfaceStatistics::computeSum(const Surface<Real> & surface){
   Real sum = 0.;
   UInt n = surface.size();
 
 #pragma omp parallel for reduction(+: sum)
   for (UInt i = 0 ; i < n*n ; ++i)
     sum += surface(i);
 
   return sum;
 }
 
 /* -------------------------------------------------------------------------- */
 std::vector<Real> SurfaceStatistics::computeMoments(Surface<Real> & surface){
 
   SurfaceComplex<Real> power = computePowerSpectrum(surface);
   Real L = surface.getL();
   UInt n = surface.size();
   Real k_factor = 2.*M_PI/L;
 
   std::vector<std::vector<Real> > moment;
   moment.resize(6);
   for (int i = 0; i < 5; i++)
     moment[i].resize(5);
 
   for (UInt i = 0; i < n/2; i++) {
     Real i0 = pow(Real(i)*k_factor,0);
     Real i2 = pow(Real(i)*k_factor,2);
     Real i4 = pow(Real(i)*k_factor,4);
     Real im0 = pow(-Real(i)*k_factor,0);
     Real im2 = pow(-Real(i)*k_factor,2);
     Real im4 = pow(-Real(i)*k_factor,4);
 
     for (UInt j = 0; j < n/2; j++) {
       Real Phi = power(i,j).real();
 
       Real Phi_10;
       Real Phi_01;
       Real Phi_11;
       if (i==0 && j==0) {
         Phi_10 = 0;
         Phi_11 = 0;
         Phi_01 = 0;
       }
       else if (i==0) {
         Phi_10 = 0;
         Phi_11 = 0;
         Phi_01 = power(i,(n-j)).real();
       }
       else if (j==0) {
         Phi_10 = power(n-i,j).real();
         Phi_11 = 0;
         Phi_01 = 0;
       }
       else {
         Phi_10 = power(n-i,j).real();
         Phi_11 = power(n-i,(n-j)).real();
         Phi_01 = power(i,(n-j)).real();
       }
       Real j0 = pow(Real(j)*k_factor,0);
       Real j2 = pow(Real(j)*k_factor,2);
       Real j4 = pow(Real(j)*k_factor,4);
       Real jm0 = pow(-Real(j)*k_factor,0);
       Real jm2 = pow(-Real(j)*k_factor,2);
       Real jm4 = pow(-Real(j)*k_factor,4);
 
       moment[0][0] += i0*j0*Phi + im0*j0*Phi_10 + i0*jm0*Phi_01 + im0*jm0*Phi_11;
       moment[0][2] += i0*j2*Phi + im0*j2*Phi_10 + i0*jm2*Phi_01 + im0*jm2*Phi_11;
       moment[2][0] += i2*j0*Phi + im2*j0*Phi_10 + i2*jm0*Phi_01 + im2*jm0*Phi_11;
       moment[2][2] += i2*j2*Phi + im2*j2*Phi_10 + i2*jm2*Phi_01 + im2*jm2*Phi_11;
       moment[0][4] += i0*j4*Phi + im0*j4*Phi_10 + i0*jm4*Phi_01 + im0*jm4*Phi_11;
       moment[4][0] += i4*j0*Phi + im4*j0*Phi_10 + i4*jm0*Phi_01 + im4*jm0*Phi_11;
     }
   }
 
   std::vector<Real> mean_moments;
   //m2
   Real mean_m2 = 0.5*(moment[0][2]+moment[2][0]);
   mean_moments.push_back(mean_m2);
   //m4
   Real mean_m4 = (3*moment[2][2]+moment[0][4]+moment[4][0])/3.;
   mean_moments.push_back(mean_m4);
   return mean_moments;
 }
 
 /* -------------------------------------------------------------------------- */
 
 Real SurfaceStatistics::computeAnalyticFractalMoment(UInt order,UInt k1, UInt k2, Real hurst,Real A, Real L){
   Real xi = k2/k1;
   Real exp = - 2.*hurst;
   Real _k1 = 2.*M_PI/L*k1;
 
   Real Cq;
   switch (order){
   case 0: Cq = 2.*M_PI;break;
   case 2: Cq = M_PI;break;
   case 4: Cq = 3./4.*M_PI;break;
   default: SURFACE_FATAL("invalid order");
   }
   exp = order - 2.*hurst;
   Real res = A * Cq * pow(_k1,exp) / exp * ( pow(xi, exp) - 1 );
 
   return res;
 }
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
 
 #endif /* __SURFACE_STATISTICS_HH__ */
diff --git a/tests/test_fftransform.py b/tests/test_fftransform.py
index 52bcd91..08e2eea 100644
--- a/tests/test_fftransform.py
+++ b/tests/test_fftransform.py
@@ -1,70 +1,70 @@
 #!/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
 
 def main():
     size = 8
     # Generate random rough surface
     SG = tm.SurfaceGeneratorFilterFFT()
     SG.getGridSize().assign(size)
     SG.getHurst().assign(0.8)
     SG.getRMS().assign(1.);
     SG.getQ0().assign(4);
     SG.getQ1().assign(4);
     SG.getQ2().assign(32);
     SG.getRandomSeed().assign(156);
     SG.Init()
     surface = SG.buildSurface()
 
     # Compute fourier transform through numpy
     numpy_spectral = np.fft.rfft2(surface)
 
     # Compute fourier tranform through tamaas
     tamaas_spectral = tm.SurfaceRealComplex(size, 1.)
-    transform = tm.FFTransformFFTWReal(size, surface, tamaas_spectral)
+    transform = tm.FFTransformFFTWReal(surface, tamaas_spectral)
     transform.forwardTransform()
 
     spectral = tm.convertGrid(tamaas_spectral)
     error_npy = np.abs(np.sum((spectral-numpy_spectral)**2))
 
     surface_copy = surface.copy()
     surface[:,:] = 0
 
     transform.backwardTransform()
     error_back = np.abs(np.sum((surface_copy-surface)**2))
 
     print error_back
     print error_npy
 
     if error_back > 1e-16 or error_npy > 1e-16:
         return 1
 
     return 0
 
 if __name__ == '__main__':
     sys.exit(main())