diff --git a/python/clusters.i b/python/clusters.i
index eb4f209..f234044 100644
--- a/python/clusters.i
+++ b/python/clusters.i
@@ -1,17 +1,11 @@
 /* -------------------------------------------------------------------------- */
 /* Wrapping contact clusters classes                                           */
 /* -------------------------------------------------------------------------- */
 %ignore operator <<;
 
-%include "contact_area.hh"
-%include "contact_cluster_collection.hh"
 %include "flood_fill.hh"
 
-%include "std_vector.i"
 %include "std_list.i"
 namespace tamaas {
-%template(VectorClusters) ::std::vector<::tamaas::ContactCluster>;
-%template(ListClusters) ::std::list<::tamaas::Cluster>;
+  %template(ListClusters) ::std::list<::tamaas::Cluster>;
 }
-%include "contact_cluster.hh"
-
diff --git a/python/surface.i b/python/surface.i
index addcdef..0d6c495 100644
--- a/python/surface.i
+++ b/python/surface.i
@@ -1,388 +1,378 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 
 %{
 #include <complex>
-#include "map_2d.hh"
 #include "surface.hh"
 #include "types.hh"
 #define NPY_NO_DEPRECATED_API NPY_1_7_API_VERSION
 #include <numpy/arrayobject.h>
 #include "surface_generator_filter_fft.hh"
 #include "smart_pointer.hh"
 #include "surface_statistics.hh"
 #include "bem_kato.hh"
 #include "bem_polonski.hh"
 #include "bem_penalty.hh"
 #include "bem_gigi.hh"
 #include "bem_gigipol.hh"
 #include "bem_uzawa.hh"
 #include "bem_fft_base.hh"
 #include "bem_grid.hh"
 #include "bem_grid_polonski.hh"
 #include "bem_grid_kato.hh"
 #include "bem_grid_teboulle.hh"
 #include "bem_grid_condat.hh"
-#include "contact_area.hh"
-#include "contact_cluster_collection.hh"
 #include "functional.hh"
 #include "meta_functional.hh"
 #include "elastic_energy_functional.hh"
 #include "complimentary_term_functional.hh"
 #include "exponential_adhesion_functional.hh"
 #include "squared_exponential_adhesion_functional.hh"
 #include "maugis_adhesion_functional.hh"
 #include "surface_timer.hh"
 #include "fftransform_fftw.hh"
 #include "flood_fill.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(unsigned_int,tamaas::UInt) // dunno why this doesn't work
 %my_pointer_class(vectorReal,std::vector<tamaas::Real>)
 
 %apply (unsigned int IN_ARRAY1[ANY]) {(tamaas::UInt n[ANY])};
 
 /* -------------------------------------------------------------------------- */
 /* Wrapping Surface class                                                     */
 /* -------------------------------------------------------------------------- */
 
 %include "tamaas.hh"
 
 %ignore tamaas::Grid::iterator::operator++;
 %ignore tamaas::Grid::begin;
 %ignore tamaas::Grid::end;
 %ignore *::operator=;
 
 %include "grid.hh"
 %include "grid_hermitian.hh"
 namespace tamaas {
   %template (Grid1dReal) Grid<Real, 1>;
   %template (Grid2dReal) Grid<Real, 2>;
   %template (Grid2dBool) Grid<bool, 2>;
   %template (Grid2dComplex) Grid<Complex, 2>;
   %template (Grid2dInt) Grid<Int, 2>;
   %template (Grid2dUInt) Grid<UInt, 2>;
   %template (GridHermitian2dReal) GridHermitian<Real, 2>;
 }
 
 %include "types.hh"
 
 %include "grid_hermitian.hh"
 
-%include "map_2d.hh"
-namespace tamaas {
-%template(Map2dComplex) Map2d<Complex>;
-%template(Map2dInt)     Map2d<int>;
-%template(Map2dReal)    Map2d<Real>;
-}
-
 %include "map_2d_square.hh"
 namespace tamaas {
 %template(Map2dSquareComplex) Map2dSquare<Complex>;
 %template(Map2dSquareInt)     Map2dSquare<int>;
 %template(Map2dSquareReal)    Map2dSquare<Real>;
 }
 
 %include "surface.hh"
 %include "surface_complex.hh"
 
 namespace tamaas {
 %template(SurfaceReal) Surface<Real>;
 %template(SurfaceInt)  Surface<int>;
 %template(SurfaceRealComplex) SurfaceComplex<Real>;
 }
 
 %inline %{
   __BEGIN_TAMAAS__
     template <int required_type>
     void sanityCheck(PyObject * input) {
       if (PyArray_TYPE((PyArrayObject*)input) != required_type)
 	throw Exception("Sanity check: python object stores wrong type");
     }
   __END_TAMAAS__
 %}
 
 %inline %{
   __BEGIN_TAMAAS__
 
   template <typename T, UInt d>
   class GridForPython : public Grid<T, d> {
     public:
     GridForPython(T * wrapped_memory, UInt n[d], Real L[d], UInt nb_components):
     Grid<T, d>() {
       this->__init__(wrapped_memory, n, nb_components);
     }
 
     void __init__(T * w, UInt n[d], UInt nb_components) {
       memcpy(this->n, n, d * sizeof(UInt));
       memset(this->L, 1, d * sizeof(Real));
       this->nb_components = nb_components;
       this->data.wrapMemory(w, this->computeSize());
     }
 
     GridForPython(PyObject * input): Grid<T, d>() {
       if (!PyArray_Check(input)) {
 	PyObject* obj_type = PyObject_Type(input);
 	std::string type = PyString_AsString(obj_type);
 	throw Exception("GridForPython: incompatible input which is not a numpy: " + type);
       }
 
       UInt _n = PyArray_NDIM((PyArrayObject*)input); // numpy number of dimensions
       npy_intp * ndims = PyArray_DIMS((PyArrayObject*)input); // numpy dimensions
       UInt grid_dims[d] = {0}; // tamaas dims
       if (_n >= d) {
 	for (UInt i = 0 ; i < d ; ++i)
 	  grid_dims[i] = ndims[i];
       } else {
 	throw Exception("GridForPython: incompatible input size for numpy array: " + _n);
       }
 
       if (_n > d) this->nb_components = ndims[d];
       else this->nb_components = 1;
 
       PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew(input);
       if (iter == NULL) throw;
       this->__init__(( T *)(iter->dataptr), grid_dims, this->nb_components);
       Py_DECREF(iter);
     }
 
     virtual ~GridForPython() {}
   };
 
   __END_TAMAAS__
 %}
 
 %inline %{
   __BEGIN_TAMAAS__
 template <typename T>
   class SurfaceForPython : public Surface<T>{
 
  public:
 
  SurfaceForPython(T * wrapped_memory,
 		   UInt n,
 		   Real L)
    : Surface<T>(0,L){
     this->__init__(wrapped_memory,n,L);
   };
 
 
   void __init__(T * wrapped_memory,
 		UInt n,
 		Real L){
     this->n[0] = n;
     this->n[1] = n;
     this->data.wrapMemory(wrapped_memory,n*n);
   }
 
   SurfaceForPython(PyObject * input, Real L = 1.)
     : Surface<T>(0,L){
 
     if (!PyArray_Check(input)){
       PyObject* obj_type = PyObject_Type(input);
       std::string type = PyString_AsString(obj_type);
       throw Exception("SurfaceForPython: incompatible input which is not a numpy: " + type);
     }
 
     UInt _n = PyArray_NDIM((PyArrayObject*)input);
     if (_n != 2) SURFACE_FATAL("SurfaceForPython: incompatible numpy dimension " << _n);
     npy_intp * ndims = PyArray_DIMS((PyArrayObject*)input);
     UInt sz = ndims[0];
     UInt sz2 = ndims[1];
     if (sz != sz2) SURFACE_FATAL("SurfaceForPython: incompatible numpy dimensions " << sz << " " << sz2);
     PyArrayIterObject *iter = (PyArrayIterObject *)PyArray_IterNew(input);
     if (iter == NULL) throw;
     this->__init__(( T *)(iter->dataptr),sz,L);
     Py_DECREF(iter);
   };
 
   ~SurfaceForPython(){
   };
 
   void resize(UInt new_size){
     if (this->size() != new_size)
       throw Exception("SurfaceForPython: cannot resize a temporary vector");
   }
  };
   __END_TAMAAS__
 %}
 
 %define %grid_typemaps(GRID_CLASS, DATA_TYPE, DATA_TYPECODE, DIM)
 %typemap(out, fragment="NumPy_Fragments") (GRID_CLASS &)
 {
   using namespace tamaas;
   UInt additional_dim = ($1->getNbComponents() == 1)? 0:1;
   UInt dim = DIM+additional_dim;
   npy_intp dims[dim];
   for (UInt i = 0 ; i < DIM ; i++) dims[i] = (npy_intp)$1->sizes()[i];
   if (additional_dim)
     dims[dim-1] = $1->getNbComponents();
 
   DATA_TYPE * data = const_cast<DATA_TYPE *>($1->getInternalData());
   PyObject* obj = PyArray_SimpleNewFromData(dim,
 					    dims,
 					    DATA_TYPECODE,
 					    &data[0]);
   PyArrayObject* array = (PyArrayObject*) obj;
   if (!array) SWIG_fail;
   $result = SWIG_Python_AppendOutput($result, obj);
 }
 
 %typemap(in) tamaas::Grid<DATA_TYPE, DIM> & {
   using namespace tamaas;
   if (!PyArray_Check($input)) {
     PyObject* objectsRepresentation = PyObject_Repr(PyObject_Type($input));
     const char* s = PyString_AsString(objectsRepresentation);
     std::string __type = s;
 
     if ("<class 'tamaas.GridForPython" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
-    else if ("<class 'tamaas.Grid" + std::string("DATA_TYPE") + "'>" == __type) {
+    else if ("<class 'tamaas.Grid" + std::string("DIM") + "d" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else throw tamaas::Exception("typemap: incompatible input which is not a proper type: " + __type + " (DATATYPE = " + std::string("DATA_TYPE") + ")");
   } else {
     $1 = new GridForPython< DATA_TYPE, DIM >($input);
     sanityCheck< DATA_TYPECODE >($input);
   }
 }
 
 %inline %{
   namespace tamaas {
   GRID_CLASS & convertGrid(GRID_CLASS & s){
     return s;
   }
   }
 %}
 %enddef
 
 
 %define %surface_typemaps(SURFACE_CLASS, DATA_TYPE, DATA_TYPECODE)
 
 %typemap(out, fragment="NumPy_Fragments") (SURFACE_CLASS &)
 {
   using namespace tamaas;
   UInt nz = $1->size();
   npy_intp dims[2] = {nz,nz};
   DATA_TYPE * data = const_cast<DATA_TYPE *>($1->getInternalData());
   PyObject* obj = PyArray_SimpleNewFromData(2, dims, DATA_TYPECODE, &data[0]);
   PyArrayObject* array = (PyArrayObject*) obj;
   if (!array) SWIG_fail;
   $result = SWIG_Python_AppendOutput($result, obj);
 }
 
 %typemap(in) tamaas::Surface< DATA_TYPE > & {
   using namespace tamaas;
   if (!PyArray_Check($input)){
     PyObject* objectsRepresentation = PyObject_Repr(PyObject_Type($input));
     const char* s = PyString_AsString(objectsRepresentation);
     std::string __type = s;
 
     if ("<class 'tamaas.SurfaceForPython" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else if ("<class 'tamaas.Surface" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else if ("<class 'tamaas.SurfaceComplex" + std::string("DATA_TYPE") + "'>" == __type) {
       SWIG_ConvertPtr($input, (void **) &$1, $1_descriptor,0);
     }
     else throw tamaas::Exception("typemap: incompatible input which is not a proper type: " + __type + " (DATATYPE = " + std::string("DATA_TYPE") + ")");
   }
   else{
     $1 = new SurfaceForPython< DATA_TYPE >($input);
     sanityCheck< DATA_TYPECODE >($input);
   }
  }
 
 %inline %{
   namespace tamaas {
   SURFACE_CLASS & convertSurface(SURFACE_CLASS & s){
     return s;
   }
   }
 %}
 
 
 %enddef
 
 %surface_typemaps(tamaas::Surface<Real>, Real,       NPY_DOUBLE)
 %surface_typemaps(tamaas::SurfaceComplex<Real>, std::complex<Real>,        NPY_COMPLEX128 )
 %surface_typemaps(tamaas::Surface<unsigned long>, unsigned long, NPY_UINT  )
 %surface_typemaps(tamaas::Surface<UInt>, UInt, NPY_UINT  )
 %surface_typemaps(tamaas::Surface<int>, int, NPY_INT  )
 
 %grid_typemaps(tamaas::Grid1dReal, Real, NPY_DOUBLE, 1)
 %grid_typemaps(tamaas::Grid2dReal, Real, NPY_DOUBLE, 2)
 %grid_typemaps(tamaas::Grid2dBool, bool, NPY_BOOL, 2)
 %grid_typemaps(tamaas::GridHermitian2dReal, tamaas::Complex, NPY_COMPLEX128, 2)
 %grid_typemaps(tamaas::Grid2dUInt, UInt, NPY_UINT, 2)
 %grid_typemaps(tamaas::SurfaceComplex<Real>, tamaas::Complex, NPY_COMPLEX128, 2)
 
 namespace tamaas {
 %template(SurfaceForPythonReal) SurfaceForPython< Real >;
 }
 
 %include "surface_timer.hh"
diff --git a/python/surface_statistics.i b/python/surface_statistics.i
index 651d118..89b3dfa 100644
--- a/python/surface_statistics.i
+++ b/python/surface_statistics.i
@@ -1,46 +1,38 @@
 /**
  *
  * @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 surface statistics
 /* -------------------------------------------------------------------------- */
 
 //%include "surface.i"
 %include "surface_statistics.hh"
 %template(computeDistribution)          tamaas::SurfaceStatistics::computeDistribution<false>;
 %template(computeDistributionWithZeros) tamaas::SurfaceStatistics::computeDistribution<true>;
 %template(vecInt) std::vector<int>;
 //%include "contact_area.hh"
 //%include "contact_cluster_collection.hh"
 
-namespace tamaas{
-%extend ContactArea
-{
-  Surface<int> & getSurface() {
-    return *$self;
-  }
-}
-}
diff --git a/src/SConscript b/src/SConscript
index f8078e6..8aa0cbe 100644
--- a/src/SConscript
+++ b/src/SConscript
@@ -1,78 +1,73 @@
 Import('main_env')
 
 import os
 
 env = main_env
 
 # Lib roughcontact
 generator_list = """
 surface_generator.cpp
 surface_generator_filter.cpp
 surface_generator_filter_fft.cpp
 isopowerlaw.cpp
 """.split()
 #env.SharedLibrary('Generator',generator_list)
 
 # Lib SURFACE
 surface_list = """
-map_2d.cpp
 map_2d_square.cpp
 surface.cpp
 surface_timer.cpp
 tamaas.cpp
 """.split()
 #env.SharedLibrary('Surface',surface_list)
 
 # Lib PERCOLATION
 percolation_list = """
-cluster_grow.cpp
-contact_area.cpp
-contact_cluster.cpp
-contact_cluster_collection.cpp
 flood_fill.cpp
 """.split()
 #env.SharedLibrary('Percolation',percolation_list)
 
 # BEM PERCOLATION
 bem_list = """
 tamaas_info.cpp
 bem_kato.cpp
 bem_polonski.cpp
 bem_gigi.cpp
 bem_gigipol.cpp
 bem_penalty.cpp
 bem_uzawa.cpp
 bem_fft_base.cpp
 functional.cpp
 meta_functional.cpp
 elastic_energy_functional.cpp
 exponential_adhesion_functional.cpp
 squared_exponential_adhesion_functional.cpp
 maugis_adhesion_functional.cpp
 complimentary_term_functional.cpp
 fftransform.cpp
 fftransform_fftw.cpp
 fft_plan_manager.cpp
 grid.cpp
 grid_hermitian.cpp
 types.cpp
 bem_grid.cpp
 bem_grid_polonski.cpp
 bem_grid_kato.cpp
 bem_grid_teboulle.cpp
 bem_grid_condat.cpp
 """.split()
 #env.SharedLibrary('BEM',bem_list)
 
 rough_contact_list = generator_list + surface_list + percolation_list + bem_list
 current_files = set(os.listdir(Dir('.').srcnode().abspath))
 unexpected_files = list(current_files - set(rough_contact_list + ['SConscript']))
 
 unexpected_headers = [ f for f in unexpected_files if os.path.splitext(f)[1] == '.hh' ]
 
 unexpected_files = [ f for f in unexpected_files if not os.path.splitext(f)[1] == '.hh' ]
 
 #print "Unexpected files : {}".format(unexpected_files)
 #print "Include path : {}".format([str(x) for x in env['CPPPATH']])
 #print unexpected_headers
 env.SharedLibrary('Tamaas',rough_contact_list)
diff --git a/src/array.hh b/src/array.hh
index 06c748a..ea7d052 100644
--- a/src/array.hh
+++ b/src/array.hh
@@ -1,137 +1,134 @@
 /**
  *
  * @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 __ARRAY__HH__
 #define __ARRAY__HH__
 /* -------------------------------------------------------------------------- */
-#include <cstring>
 #include "tamaas.hh"
-
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 template <typename T>
 class Array {
 public:
   /// Copy constructor (deep)
   Array(const Array & v):
     _data(NULL),
     _size(0),
     wrapped(false) {
     this->resize(v._size);
-    memcpy(this->_data,v.getPtr(),this->_size*sizeof(T));
+    std::copy(v.getPtr(), v.getPtr()+v._size, this->_data);
   }
 
   Array():
     _data(NULL),
     _size(0),
     wrapped(false) {}
 
   Array(UInt size):
     _data(NULL),
     _size(0),
     wrapped(false) {
     this->resize(size);
     this->wrapped = false;
   }
 
   Array(T * data, UInt size):
     _data(NULL),
     _size(0),
     wrapped(false) {
     this->wrapMemory(data,size);
   }
 
   virtual ~Array(){
     if (wrapped == false && _data != NULL){
       delete [] _data;
     }
   }
 
   Array & operator=(const Array & v){
     this->wrapped = false;
     if (v.size() != this->size()) this->resize(v.size());
     memcpy(this->_data,v.getPtr(),this->_size*sizeof(T));
     return *this;
   }
 
   void wrapMemory(T * data, UInt size) {
     this->_data = data;
     this->_size = size;
     this->wrapped = true;
   }
 
   void assign(UInt size, const T & value){
     this->resize(size);
     for (UInt i = 0; i < size; i++) {
       _data[i] = value;
     }
   }
 
   const T * getPtr() const{ return _data; }
   T * getPtr() { return _data; }
 
   void setPtr(T * new_ptr) { _data = new_ptr;}
 
   void resize(UInt size){
     if (size == 0 && wrapped == false) {
       if (_data != NULL) {
 	delete [] this->_data;
 	this->_data = NULL;
       }
       this->_size = 0;
       return;
     }
     if (size == this->_size) return;
     if (wrapped == true) TAMAAS_EXCEPTION("cannot resize wrapped surface");
     if (_data != NULL) delete [] _data;
     this->_data = new T[size];
     this->_size = size;
   }
 
   inline T & operator[](UInt i){
     TAMAAS_ASSERT(i < _size, "memory overflow");
     return _data[i];
   }
 
   inline const T & operator[](UInt i)const{
     TAMAAS_ASSERT(i < _size, "memory overflow");
     return _data[i];
   }
 
   inline UInt size() const{return _size;}
 
 private:
   T * _data;
   UInt _size;
 
   bool wrapped;
 };
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 #endif /* __ARRAY__HH__ */
diff --git a/src/cluster_grow.cpp b/src/cluster_grow.cpp
deleted file mode 100644
index 0be39b2..0000000
--- a/src/cluster_grow.cpp
+++ /dev/null
@@ -1,107 +0,0 @@
-/**
- *
- * @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 "cluster_grow.hh"
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-ClusterGrow::ClusterGrow(Map2d<int> & area,
-			 int initial_state,
-			 int to_check_state,
-			 int contact_state,
-			 int no_contact_state):
-  current_state   (initial_state),
-  to_check_state  (to_check_state),
-  contact_state   (contact_state),
-  no_contact_state(no_contact_state),
-  area(area){
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-ClusterGrow::~ClusterGrow(){
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-
-void ClusterGrow::addIndicesToExplore(std::vector<std::pair<UInt,UInt> > &
-				      indices_to_add){
-
-  auto it  = indices_to_add.begin();
-  auto end = indices_to_add.end();
-
-  for( ; it != end ; ++it){
-    int & val = area(*it);
-
-    if (val == contact_state) {
-
-      index_to_explore.push(*it);
-
-      val = to_check_state;
-    }
-  }
-}
-/* -------------------------------------------------------------------------- */
-
-template <bool periodic>
-void ClusterGrow::grow(UInt i, UInt j){
-
-  if (area(i,j) != contact_state) return;
-
-  ++current_state;
-
-  area(i,j) = current_state;
-
-  std::vector<std::pair<UInt,UInt> > v = area.getNeighborIndexes<periodic>(i,j);
-  this->addIndicesToExplore(v);
-
-  while (!index_to_explore.empty()){
-
-    std::pair<UInt,UInt> ind = index_to_explore.front();
-
-    int & val = area(ind);
-
-    if (val == contact_state || val == to_check_state){
-      val = current_state;
-      std::vector<std::pair<UInt,UInt> > v = area.getNeighborIndexes<periodic>(ind);
-      this->addIndicesToExplore(v);
-    }
-
-    index_to_explore.pop();
-
-  }
-}
-
-/* -------------------------------------------------------------------------- */
-
-template void ClusterGrow::grow<true>(UInt i, UInt j);
-template void ClusterGrow::grow<false>(UInt i, UInt j);
-
-__END_TAMAAS__
diff --git a/src/cluster_grow.hh b/src/cluster_grow.hh
deleted file mode 100644
index 31b6d7c..0000000
--- a/src/cluster_grow.hh
+++ /dev/null
@@ -1,105 +0,0 @@
-/**
- *
- * @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 __CLUSTER_GROW_HH__
-#define __CLUSTER_GROW_HH__
-/* -------------------------------------------------------------------------- */
-#include "contact_area.hh"
-/* -------------------------------------------------------------------------- */
-#include <queue>
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-
-class ClusterGrow {
-
-  /* ------------------------------------------------------------------------ */
-  /* Typedefs                                                                 */
-  /* ------------------------------------------------------------------------ */
-
-public:
-
-  /* ------------------------------------------------------------------------ */
-  /* Constructors/Destructors                                                 */
-  /* ------------------------------------------------------------------------ */
-public:
-
-  ClusterGrow(Map2d<int> & area,
-	      int initial_state,
-	      int to_check_state = -1,
-	      int contact_state = 1,
-	      int no_contact_state = 0);
-  virtual ~ClusterGrow();
-
-  /* ------------------------------------------------------------------------ */
-  /* Methods                                                                  */
-  /* ------------------------------------------------------------------------ */
-public:
-
-  //! start grow of a cluster from a given point
-  template <bool periodic>
-  void grow(UInt i, UInt j);
-
-private:
-
-  void addIndicesToExplore(std::vector<std::pair<UInt,UInt> > & indices_to_add);
-
-  /* ------------------------------------------------------------------------ */
-  /* Accessors                                                                */
-  /* ------------------------------------------------------------------------ */
-
-public:
-
-  UInt getNbClusters(){
-    UInt tmp = current_state - 2;
-    return tmp;
-  };
-
-
-  /* ------------------------------------------------------------------------ */
-  /* Class Members                                                            */
-  /* ------------------------------------------------------------------------ */
-
-private:
-
-  int current_state;
-
-  const int to_check_state  ;
-  const int contact_state   ;
-  const int no_contact_state;
-
-  Map2d<int> & area;
-
-  std::queue<std::pair<UInt,UInt> > index_to_explore;
-};
-
-/* -------------------------------------------------------------------------- */
-
-__END_TAMAAS__
-
-
-#endif /* __CLUSTER_GROW_HH__ */
diff --git a/src/contact_area.cpp b/src/contact_area.cpp
deleted file mode 100644
index fbdd3bb..0000000
--- a/src/contact_area.cpp
+++ /dev/null
@@ -1,146 +0,0 @@
-/**
- *
- * @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"contact_area.hh"
-#include <iostream>
-#include "cluster_grow.hh"
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-ContactArea::ContactArea(UInt n, Real L): Surface<int>(n,L) {
-  //num_cluster = 0;
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <>
-ContactArea::ContactArea(const Map2dSquare<Real> & m) :
-  Surface<int>(m.size(),m.getL()){
-
-  UInt n = this->size();
-
-  for (UInt i = 0 ; i < n*n ; ++i)
-    this->at(i) = (m(i) > 0.);
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-
-ContactArea::~ContactArea() {}
-
-/* -------------------------------------------------------------------------- */
-
-void ContactArea::detectContactClusters() {
-
-  // starting state was push as 2 by vlad
-  ClusterGrow clg(*this,2);
-  UInt N = this->size();
-
-  for (UInt i = 0; i < N; ++i) {
-      for (UInt j = 0; j < N; ++j) {
-	clg.grow<true>(i,j);
-      }
-  }
-
-  for (UInt i = 0; i < N; ++i) {
-      for (UInt j = 0; j < N; ++j) {
-	if (this->at(i,j) == -1)
-	  SURFACE_FATAL("problem " << i << " " << j << " is -1");
-      }
-  }
-
-  for (UInt i = 0; i < clg.getNbClusters(); ++i)
-    colors.push_back(i+3);
-}
-/* -------------------------------------------------------------------------- */
-
-
-void ContactArea::replaceColors(int from, int to) {
-  // Repalce color "from" by the color "to"
-  UInt n = this->size();
-  for (UInt i = 0; i < n; i++)
-    for (UInt j = 0; j < n; j++)
-      if (this->at(i,j) == from)
-	this->at(i,j) = to;
-}
-/* -------------------------------------------------------------------------- */
-
-void ContactArea::randomizeClusterID(UInt seed) {
-  std::cout << "/ Randomize colors of contact clusters." << std::endl;
-  UInt shift=1000, tryid;
-  bool found;
-  std::vector<UInt> cid;
-  UInt num_cluster = colors.size();
-  cid.resize(num_cluster);
-  while (shift < num_cluster)
-    shift*=10;
-  srand(seed);
-  for (int i = 0; i < int(colors.size()); i++) {
-    tryid=-1;
-    do {
-      tryid = shift+int((num_cluster+10)*random()/Real(RAND_MAX));
-      found = true;
-      for (int j = 0; j < i; j++)
-	if (cid[j] == tryid)
-          found = false;
-    }
-    while (!found);
-    cid[i]=tryid;
-    replaceColors(colors[i],tryid);
-    colors[i]=tryid;
-  }
-  // Shift all colors back to 1
-  for (int i = 0; i < int(colors.size()); i++) {
-    replaceColors(colors[i],colors[i]-shift+1);
-    colors[i]-=shift-1;
-  }
-}
-
-/* -------------------------------------------------------------------------- */
-
-std::vector<ContactCluster> ContactArea::getContactClusters() {
-
-  std::vector<ContactCluster> clusters;
-  UInt nb_colors = this->colors.size();
-
-  for (UInt i = 0; i < nb_colors; i++) {
-    //    std::cout << "Extracting contact cluster " << i << std::endl;
-    ContactCluster c(colors[i], *this);
-    //std::cout << "Compute statistics of contact cluster " << i << std::endl;
-    c.computeStatistics();
-    clusters.push_back(c);
-  }
-
-  return clusters;
-}
-
-/* -------------------------------------------------------------------------- */
-
-
-
-__END_TAMAAS__
diff --git a/src/contact_area.hh b/src/contact_area.hh
deleted file mode 100644
index 6528c4d..0000000
--- a/src/contact_area.hh
+++ /dev/null
@@ -1,70 +0,0 @@
-//==================================================================//
-//  Percolation detection and characterization of contact clusters  //
-//  To be integrated in FFT-BEM code                                //
-//                                               		    //
-//  V.A. Yastrebov, 2014                         		    //
-//  Centre des Materiaux, CNRS, MINES ParisTech  		    //
-//								    //
-//  G. Anciaux, 2014						    //
-//  EPFL, ENAC, IIC, LSMS					    //
-//==================================================================//
-
-#ifndef __CONTACT_AREA__
-#define __CONTACT_AREA__
-/* -------------------------------------------------------------------------- */
-#include "map_2d.hh"
-#include <vector>
-#include "contact_cluster.hh"
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-
-class ContactArea : public Surface<int> {
-
-  /* ------------------------------------------------------------------------ */
-  /* Constructors/Destructors                                                 */
-  /* ------------------------------------------------------------------------ */
-
-public:
-
-  ContactArea(UInt n, Real L);
-
-  template <typename T>
-  ContactArea(const Map2dSquare<T> & n);
- ~ContactArea();
-
-  /* ------------------------------------------------------------------------ */
-  /* Methods                                                                  */
-  /* ------------------------------------------------------------------------ */
-
-public:
-
-  //! TODO
-  void detectContactClusters();
-  //! TODO
-  void replaceColors(int from, int to);
-  //! TODO
-  void randomizeClusterID(UInt seed);
-  //! generate contact clusters
-  std::vector<::tamaas::ContactCluster> getContactClusters();
-
-  /* ------------------------------------------------------------------------ */
-  /* Accessors                                                                */
-  /* ------------------------------------------------------------------------ */
-
-  UInt getNumClusters(){return colors.size();};
-
-
-  /* ------------------------------------------------------------------------ */
-  /* Class Members                                                            */
-  /* ------------------------------------------------------------------------ */
-
-private:
-  std::vector<int> colors;
-
-};
-
-__END_TAMAAS__
-
-#endif
diff --git a/src/contact_cluster.cpp b/src/contact_cluster.cpp
deleted file mode 100644
index 85bb6c2..0000000
--- a/src/contact_cluster.cpp
+++ /dev/null
@@ -1,405 +0,0 @@
-//==================================================================//
-//  Percolation detection and characterization of contact clusters  //
-//  To be integrated in FFT-BEM code                                //
-//                                               		    //
-//  V.A. Yastrebov, 2014                         		    //
-//  Centre des Materiaux, CNRS, MINES ParisTech  		    //
-//								    //
-//  G. Anciaux, 2014						    //
-//  EPFL, ENAC, IIC, LSMS					    //
-//==================================================================//
-
-
-#include "contact_cluster.hh"
-#include <iostream>
-#include <cmath>
-#include "contact_area.hh"
-#include "cluster_grow.hh"
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-void ContactCluster::computeBoundingBox(ContactArea & contact_area) {
-
-  int N = contact_area.size();
-
-  int i0 = -1;
-  int j0 = -1;
-  int i1 = -1;
-  int j1 = -1;
-
-  for (int i=0; i<N; i++) {
-    for (int j=0; j<N; j++) {
-      if (contact_area(i,j) == this->id && i0 < 0) 	    i0 = i;
-      if (contact_area(j,i) == this->id && j0 < 0)          j0 = i;
-      if (contact_area(N-i-1,N-j-1) == this->id && i1 < 0)  i1 = N-i-1;
-      if (contact_area(N-j-1,N-i-1) == this->id && j1 < 0)  j1 = N-i-1;
-
-      // if all extremum are found then break
-      if (i0 >= 0 && j0 >= 0 && i1 >= 0 && j1 >= 0) break;
-    }
-    // if all extremum are found then break
-    if (i0 >= 0 && j0 >= 0 && i1 >= 0 && j1 >= 0)   break;
-  }
-
-  if (i0 < 0 || j0 < 0) SURFACE_FATAL("bounds not found " << id << " " << contact_area(0,0));
-  if (i1 < 0 || j1 < 0) SURFACE_FATAL("bounds not found");
-
-  lb_corner[0] = i0;
-  lb_corner[1] = j0;
-
-  rt_corner[0] = i1;
-  rt_corner[1] = j1;
-
-  bbox_size[0] = i1-i0+1;
-  bbox_size[1] = j1-j0+1;
-
-}
-
-
-/* -------------------------------------------------------------------------- */
-
-void ContactCluster::checkLines(UInt dir,ContactArea & contact_area){
-
-  UInt N  = contact_area.size();
-  UInt i0 = lb_corner[dir];
-
-  start[dir] = -1;
-  stop[dir] = -1;
-
-  for (UInt i = 0; i < N; i++) {
-
-    bool fail = false;
-    for (UInt j = 0; j < N; j++){
-      int val = -1;
-      if (dir == 0) val = contact_area(i,j);
-      if (dir == 1) val = contact_area(j,i);
-
-      if (val == id) {
-	fail = true;
-	break;
-      }
-    }
-
-    if (i0 == 0 && !fail && start[dir] < 0) {
-      start[dir] = i;
-      if (stop[dir] >= 0) SURFACE_FATAL("Error occured at "<<__LINE__);
-    }
-    if (i0 == 0 && fail && start[dir] >= 0 && stop[dir] < 0)
-      stop[dir] = i;
-  }
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-
-void ContactCluster::minimizeBoundingBox(ContactArea & contact_area) {
-
-  this->split_cluster = false;
-  UInt N = contact_area.size();
-  if (bbox_size[0] != N && bbox_size[1] != N) return;
-
-  // Minimize bounding box
-  this->split_cluster = true;
-
-  // Check if there are vertical lines that do not contact "id"
-  this->checkLines(0,contact_area);
-  // Check if there are horizontal lines that do not contact "id"
-  this->checkLines(1,contact_area);
-
-  if (start[0] >= 0) bbox_size[0] -= (stop[0]-start[0])-1;
-  if (start[1] >= 0) bbox_size[1] -= (stop[1]-start[1])-1;
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-void ContactCluster::extractCluster(ContactArea & contact_area) {
-
-  UInt Nx = this->size(0);
-  UInt Ny = this->size(1);
-  UInt i0 = lb_corner[0];
-  UInt j0 = lb_corner[1];
-  UInt N = contact_area.size();
-
-  mass_center[0] = 0;
-  mass_center[1] = 0;
-
-  for (UInt i = 0; i < Nx; i++) {
-    for (UInt j = 0; j < Ny; j++) {
-      if (i+i0 >= N || j+j0 >= N)
-	SURFACE_FATAL("    i0=" << i0
-		      << ",j0=" << j0
-		      << ", i=" <<  i
-		      << ", j=" <<  j);
-
-      int val = contact_area(i+i0,j+j0);
-      if (val == id) {
-	this->at(i,j) = 1;
-	mass_center[0] += i+i0;
-	mass_center[1] += j+j0;
-      }
-    }
-  }
-}
-
-/* -------------------------------------------------------------------------- */
-
-void ContactCluster::extractClusterPeriodic(ContactArea & contact_area) {
-
-  // if the cluster is split by periodicity
-  mass_center[0] = 0;
-  mass_center[1] = 0;
-  UInt Nx = this->size(0);
-  UInt Ny = this->size(1);
-  UInt N = contact_area.size();
-  UInt i0 = lb_corner[0];
-  UInt j0 = lb_corner[1];
-
-  UInt Nxmax = start[0]<0?Nx:N;
-  UInt Nymax = start[1]<0?Ny:N;
-
-  for (UInt i = 0; i < Nxmax; i++) {
-    for (UInt j = 0; j < Nymax; j++) {
-      if (i+i0>=0 && i+i0<N &&
-	  j+j0>=0 && j+j0<N) {
-
-	int val = contact_area(i+i0,j+j0);
-	if (val == id) {
-	  if (start[0] >= 0 && start[1] < 0) {
-	    mass_center[1] += j+j0;
-	    if ((int)i < start[0]) {
-	      this->at(i,j) = 1;
-	      mass_center[0]+=i;
-	    }
-	    else {
-	      this->at(i-(stop[0]-start[0])+1,j) = 1;
-	      mass_center[0] -= (N-i);
-	    }
-	  }
-	  else if (start[1] >= 0 && start[0] < 0) {
-	    mass_center[0]+=i+i0;
-	    if ((int)j < start[1]) {
-	      this->at(i,j) = 1;
-	      mass_center[1]+=j;
-	    }
-	    else {
-	      this->at(i,j-(stop[1]-start[1])+1) = 1;
-	      mass_center[1] -= (N-j);
-	    }
-	  }
-	  else if (start[1] >= 0 && start[0] >= 0) {
-	    if ((int)j < start[1] && (int)i < start[0]) {
-	      this->at(i,j) = 1;
-	      mass_center[0] += i;
-	      mass_center[1] += j;
-	    }
-	    else if ((int)j < start[1] && (int)i >= start[0]) {
-	      this->at(i-(stop[0]-start[0])+1,j) = 1;
-	      mass_center[1]+=j;
-	      mass_center[0]-=(N-i);
-	    }
-	    else if ((int)j >= start[1] && (int)i >= start[0]) {
-	      this->at(i-(stop[0]-start[0])+1,j-(stop[1]-start[1])+1) = 1;
-	      mass_center[0] -= (N-i);
-	      mass_center[1] -= (N-j);
-	    }
-	    else if ((int)j >= start[1] && (int)i < start[0]) {
-	      this->at(i,j-(stop[1]-start[1])+1) = 1;
-	      mass_center[0] += i;
-	      mass_center[1] -= (N-j);
-	    }
-	  }
-	  else if (start[0] < 0 && start[1] < 0) {
-	    this->at(i,j) = 1;
-	    mass_center[0] += i+i0;
-	    mass_center[1]+=j+j0;
-	  }
-	}
-      }
-      else {
-	std::cout<<"Error i0="<<i0<<",j0="<<j0
-		 <<",i="<<i<<",j="<<j
-		 <<",Nxmax="<<Nxmax << ",Nymax" << Nymax
-		 <<std::endl;
-	abort();
-      }
-    }
-  }
-}
-
-/* -------------------------------------------------------------------------- */
-
-
-ContactCluster::ContactCluster(int color, ContactArea & contact_area)
-  : Map2d(0,1.)  {
-
-  concave = false;
-  id = 0;
-  lb_corner[0] = 0.;
-  lb_corner[1] = 0.;
-  rt_corner[0] = 0.;
-  rt_corner[1] = 0.;
-  start[0] = 0.;
-  stop[1]  = 0.;
-  bbox_size[0] = 0.;
-  bbox_size[1] = 0.;
-  A = 0.;
-  P = 0.;
-  IP = 0.;
-  mass_center[0] = 0;
-  mass_center[1] = 0;
-  Ps = 0.;
-  Cd = 0.;
-  //RI = 0.;
-  comp = 0.;
-  split_cluster = 0.;
-
-  this->id = color;
-
-  this->computeBoundingBox(contact_area);
-
-  if (bbox_size[0] < 0 || bbox_size[1] < 0) {
-    SURFACE_FATAL("Something wrong with the detection of the bounding box of a contact cluster");
-  }
-
-  this->minimizeBoundingBox(contact_area);
-  this->setGridSize(bbox_size);
-
-  if (!split_cluster) this->extractCluster(contact_area);
-  else this->extractClusterPeriodic(contact_area);
-}
-
-/* -------------------------------------------------------------------------- */
-
-void ContactCluster::computeStatistics() {
-
-  bool start = false;
-  //  UInt start_i, start_j;
-  A=0;
-  P=0;
-  IP=0;
-
-  //  UInt Nmax = bbox_size[0] > bbox_size[1] ? bbox_size[0]:bbox_size[1];
-  UInt intersection[2] = {0,0};
-
-  for (UInt i = 0; i < bbox_size[0]; i++) {
-    intersection[0] = 0;
-    intersection[1] = 0;
-    for (UInt j = 0; j < bbox_size[1]; j++) {
-      // Area
-      A += this->at(i,j);
-
-      // Perimeter at edges
-      if ( (i==0 || i== bbox_size[0] - 1) && this->at(i,j) == 1) {
-	if (!start) {
-	  start=true;
-	  //      start_i=i;
-          //start_j=j;
-        }
-        P++;
-      }
-      if ( (j==0 || j== bbox_size[0] - 1) &&
-	   this->at(i,j) == 1) {
-	if (!start) {
-	  start=true;
-          //start_i=i;
-          //start_j=j;
-        }
-        P++;
-      }
-      // Perimeter along lines
-      if (i>0 && this->at(i,j) != this->at(i-1,j)) {
-	if (!start) {
-	  start=true;
-          //start_i=i;
-          //start_j=j;
-        }
-        P++;
-      }
-      // Perimeter along columns
-      if (j>0 && this->at(i,j) != this->at(i,j-1)) {
-        P++;
-      }
-      // Interpixel perimeter along lines
-      if (i>0 && this->at(i,j) == this->at(i-1,j) && this->at(i,j) > 0)
-	IP++;
-      // Interpixel perimeter along columns
-      if (j>0 && this->at(i,j) == this->at(i,j-1) && this->at(i,j) > 0)
-	IP++;
-    }
-  }
-
-  for (UInt k = 0; k < 2; ++k) {
-    mass_center[k] /= Real(A);
-    //if (mass_center[k] < 0) {
-    //      mass_center[k] += N;
-    //}
-  }
-
-  this->comp = P/sqrt(Real(A));
-  Cd = (A-P/4.)/(A-sqrt(1.*A));
-  if (intersection[0] > 2 || intersection[1] > 2)
-    concave=true;
-
-  this->findHoles();
-}
-
-/* -------------------------------------------------------------------------- */
-
-void ContactCluster::findHoles() {
-
-  UInt Nx = this->size(0);
-  UInt Ny = this->size(1);
-
-  ClusterGrow clg(*this,id,-1,0,id);
-
-  for (UInt i0 = 0; i0 < Nx; ++i0) {
-    for (UInt j0 = 0; j0 < Ny; ++j0) {
-
-      // check if the point touch the border and is not in contact
-      if ( (i0 == 0    || j0 == 0 ||
-	    i0 == Nx-1 || j0 == Ny-1) &&
-	   this->at(i0,j0) == 0) {
-
-	clg.grow<true>(i0,j0);
-      }
-    }
-  }
-  // find the holes
-  std::map<UInt,Real> holes_areas;
-  for (UInt i=0; i<Nx; i++)
-    for (UInt j=0; j<Ny; j++)
-      if (this->at(i,j) > id){
-	if (holes_areas.count(this->at(i,j)) == 0)
-	  holes_areas[this->at(i,j)] = 0;
-	holes_areas[this->at(i,j)] += 1;
-      }
-
-  //  this->RI = (Nx*Ny-Ai-A)/Real(A);
-}
-
-
-/* -------------------------------------------------------------------------- */
-
-
-void ContactCluster::printself(std::ostream & stream, int indent) const {
-
-  stream << mass_center[0] << " " << mass_center[1] << " "
-	 << bbox_size[0]   << " " << bbox_size[1] << " "
-	 << sqrt(A) << " "
-	 << P << " "
-	 << comp << " "
-	 << IP << " "
-	 << Cd << " "
-	 << this->getNbHoles() << " "
-    //	 << RI << " "
-	 << id;
-}
-
-/* -------------------------------------------------------------------------- */
-
-
-
-__END_TAMAAS__
diff --git a/src/contact_cluster.hh b/src/contact_cluster.hh
deleted file mode 100644
index eaa1d79..0000000
--- a/src/contact_cluster.hh
+++ /dev/null
@@ -1,167 +0,0 @@
-//==================================================================//
-//  Percolation detection and characterization of contact clusters  //
-//  To be integrated in FFT-BEM code                                //
-//                                               		    //
-//  V.A. Yastrebov, 2014                         		    //
-//  Centre des Materiaux, CNRS, MINES ParisTech  		    //
-//								    //
-//  G. Anciaux, 2014						    //
-//  EPFL, ENAC, IIC, LSMS					    //
-//==================================================================//
-
-#ifndef __CONTACT_CLUSTER__
-#define __CONTACT_CLUSTER__
-/* -------------------------------------------------------------------------- */
-#include <map>
-#include "surface.hh"
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-class ContactArea;
-class ContactClusterCollection;
-/* -------------------------------------------------------------------------- */
-
-
-/* Class Contact Cluster class
-
-   based on the publication:
-
-   [1] Ernesto Bribiesca,
-   "An easy measure of compactness for 2D and 3D shapes",
-   Pattern Recognition 41 (2008) 543--554.
-   doi:10.1016/j.patcog.2007.06.029
-
-*/
-
-class ContactCluster : public Map2d<int> {
-
-
-  /* ------------------------------------------------------------------------ */
-  /* Typedefs                                                                 */
-  /* ------------------------------------------------------------------------ */
-
-  friend class ContactClusterCollection;
-  /* ------------------------------------------------------------------------ */
-  /* Constructors/Destructors                                                 */
-  /* ------------------------------------------------------------------------ */
-
-public:
-
-  ContactCluster() : Map2d(0,1.) {}
- ~ContactCluster() {}
-  ContactCluster(int color, ContactArea & contact_area);
-
-  /* ------------------------------------------------------------------------ */
-  /* Methods                                                                  */
-  /* ------------------------------------------------------------------------ */
-
-public:
-
-  void computeBoundingBox(ContactArea & contact_area);
-  void minimizeBoundingBox(ContactArea & contact_area);
-  void extractCluster(ContactArea & contact_area);
-  void extractClusterPeriodic(ContactArea & contact_area);
-  void computeStatistics();
-  void findHoles();
-
-  void checkLines(UInt dir,ContactArea & contact_area);
-
-
-
-  //! function to print the contain of the class
-  virtual void printself(std::ostream & stream, int indent = 0) const;
-
-  /* ------------------------------------------------------------------------ */
-  /* Accessors                                                                */
-  /* ------------------------------------------------------------------------ */
-
-public:
-
-  std::vector<Real> getHolesAreas() const{
-    std::vector<Real> areas;
-    for (auto & x: holes_areas){
-      areas.push_back(x.second);
-    }
-    return areas;
-  }
-  UInt getA(){return A;};
-  UInt getP(){return P;};
-  UInt getNbHoles() const {return holes_areas.size();};
-
-  /* ------------------------------------------------------------------------ */
-  /* Class Members                                                            */
-  /* ------------------------------------------------------------------------ */
-
-private:
-
-  /* If the shape is for sure concave then concave=true,
-     however it may be concave even if concave=false.
-  */
-  //! concavity flag
-  bool concave;
-
-  //! color id of the cluster
-  int id;
-
-  //! coordinates of the left bottom corner of the bounding box (i0,j0)
-  UInt lb_corner[2];
-  //! coordinates of the right top corner of the bounding box (i1,j1)
-  UInt rt_corner[2];
-
-  //! from clusters connected through periodicity start stop is identifying inner part
-  int start[2];
-  int stop[2];
-
-
-  //! size of the bounding box (Nx,Ny)
-  UInt bbox_size[2];
-
-  //! contact area (A)
-  UInt A;
-
-  //! contact perimeter computed discretely (P)
-  UInt P;
-
-  //! interpixel perimeter (IP), the same as contact perimeter in notions of [1]
-  int IP;
-
-  //! the holes areas
-  std::map<UInt,Real> holes_areas;
-
-  //! center of mass (ic,jc)
-  Real mass_center[2];
-
-  //! corrected contact perimeter (Ps)
-  Real Ps;
-
-  //! discrete compactness (Cd=(n-P/4)/(n-sqrt(n)) ) from [1]
-  Real Cd;
-
-  /* Ratio of the internal (closed gap) area to the contact area
-     $$RI = (A0-Ai-A)/Real(A)$$,
-     where Ai is the area within the bounding box that surrounds contact cluster
-  */
-  //  Real RI;
-
-  //! compactness (comp = Ps/sqrt(A))
-  Real comp;
-
-
-  //! TODO
-  bool split_cluster;
-
-};
-
-/* -------------------------------------------------------------------------- */
-
-inline std::ostream & operator <<(std::ostream & stream, const ContactCluster & _this)
-{
-  _this.printself(stream);
-  return stream;
-}
-
-
-__END_TAMAAS__
-
-#endif
diff --git a/src/contact_cluster_collection.cpp b/src/contact_cluster_collection.cpp
deleted file mode 100644
index b66aa92..0000000
--- a/src/contact_cluster_collection.cpp
+++ /dev/null
@@ -1,234 +0,0 @@
-/**
- *
- * @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 "contact_cluster_collection.hh"
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-
-ContactClusterCollection::ContactClusterCollection(ContactArea & area):
-  area(area){
-
-  clusters = area.getContactClusters();
-
-  TotalArea = 0;
-  TotalPerimeter = 0;
-  NumberOfClustersWithHoles = 0;
-
-  for (UInt i = 0; i < clusters.size(); ++i) {
-    TotalArea += clusters[i].getA();
-    TotalPerimeter += clusters[i].getP();
-    NumberOfClustersWithHoles += clusters[i].getNbHoles();
-  }
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-ContactClusterCollection::~ContactClusterCollection(){
-
-}
-
-
-/* -------------------------------------------------------------------------- */
-
-void ContactClusterCollection::dumpAllClusters(const std::string & prefix){
-
-  for (UInt i = 0; i < clusters.size(); i++) {
-    std::stringstream sstr;
-    sstr << prefix << "-" << i << ".csv";
-    clusters[i].dumpToTextFile(sstr.str());
-  }
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-
-
-void ContactClusterCollection::dumpClusterInfoToFile(const std::string & filename){
-
-  std::ofstream clusterData(filename.c_str());
-  // Output info on every cluster separately
-  clusterData << "# N=" << area.size() << std::endl;
-  clusterData << "# Number of clusters=" << clusters.size() << std::endl;
-  clusterData << "# [1]Num [2]ic [3]jc [4]Nx [5]Ny [6]sqrt(A) [7]P [8]comp [9]IP [10]Cd [11]holes [12]RI [13]color [14]TotalNumberOfClusters [15]TotalArea [16]TotalPerimeter [17]NumberOfClustersWithHoles" << std::endl;
-  for (UInt i = 0; i < clusters.size(); i++) {
-    clusterData << i+1 << " " << clusters[i] << " "
-		<< clusters.size() << " "
-                << TotalArea << " "
-		<< TotalPerimeter << " "
-		<< NumberOfClustersWithHoles << std::endl;
-  }
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-Real ContactClusterCollection::getTotalArea(){
-  return TotalArea;
-}
-
-/* -------------------------------------------------------------------------- */
-
-
-Real ContactClusterCollection::getTotalPerimeter(){
-  return TotalPerimeter;
-}
-
-/* -------------------------------------------------------------------------- */
-
-Real ContactClusterCollection::getNbClustersWithHole(){
-  return NumberOfClustersWithHoles;
-}
-
-/* -------------------------------------------------------------------------- */
-
-Real ContactClusterCollection::getNbClusters(){
-  return clusters.size();
-}
-
-/* -------------------------------------------------------------------------- */
-
-std::vector<Real> ContactClusterCollection::getAreas(){
-
-  std::vector<Real> res;
-  for (UInt i = 0; i < clusters.size(); i++) {
-    //    std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl;
-    res.push_back(clusters[i].A);
-  }
-  return res;
-}
-
-/* -------------------------------------------------------------------------- */
-
-std::vector<Real> ContactClusterCollection::getMassCenters(UInt i){
-
-  std::vector<Real> res;
-  for (UInt i = 0; i < clusters.size(); i++) {
-    std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl;
-    res.push_back(clusters[i].mass_center[i]);
-  }
-  return res;
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-std::vector<Real> ContactClusterCollection::getPerimeters(){
-
-  std::vector<Real> res;
-  for (UInt i = 0; i < clusters.size(); i++) {
-    std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl;
-    res.push_back(clusters[i].P);
-  }
-  return res;
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-std::vector<Real> ContactClusterCollection::getCompactnesses(){
-
-  std::vector<Real> res;
-  for (UInt i = 0; i < clusters.size(); i++) {
-    std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl;
-    res.push_back(clusters[i].comp);
-  }
-  return res;
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-std::vector<Real> ContactClusterCollection::getInterPixelPerimeters(){
-
-  std::vector<Real> res;
-  for (UInt i = 0; i < clusters.size(); i++) {
-    std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl;
-    res.push_back(clusters[i].IP);
-  }
-  return res;
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-std::vector<Real> ContactClusterCollection::getDiscretePerimeters(){
-
-  std::vector<Real> res;
-  for (UInt i = 0; i < clusters.size(); i++) {
-    std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl;
-    res.push_back(clusters[i].Cd);
-  }
-  return res;
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-std::vector<Real> ContactClusterCollection::getNbHoles(){
-
-  std::vector<Real> res;
-  for (UInt i = 0; i < clusters.size(); i++) {
-    std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl;
-    res.push_back(clusters[i].getNbHoles());
-  }
-  return res;
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-// std::vector<Real> ContactClusterCollection::getRIs(){
-//   std::vector<Real> res;
-//   for (UInt i = 0; i < clusters.size(); i++) {
-//     std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl;
-//     res.push_back(clusters[i].RI);
-//   }
-//   return res;
-
-// }
-
-/* -------------------------------------------------------------------------- */
-
-std::vector<Real> ContactClusterCollection::getColors(){
-  std::vector<Real> res;
-  for (UInt i = 0; i < clusters.size(); i++) {
-    std::cerr << "AAAAAAAA " <<clusters[i].mass_center[i] << std::endl;
-    res.push_back(clusters[i].id);
-  }
-  return res;
-
-}
-
-/* -------------------------------------------------------------------------- */
-
-
-
-__END_TAMAAS__
diff --git a/src/contact_cluster_collection.hh b/src/contact_cluster_collection.hh
deleted file mode 100644
index 54d5592..0000000
--- a/src/contact_cluster_collection.hh
+++ /dev/null
@@ -1,97 +0,0 @@
-/**
- *
- * @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 __CONTACT_CLUSTER_COLLECTION_HH__
-#define __CONTACT_CLUSTER_COLLECTION_HH__
-/* -------------------------------------------------------------------------- */
-#include "contact_area.hh"
-#include "contact_cluster.hh"
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-
-class ContactClusterCollection {
-  /* ------------------------------------------------------------------------ */
-  /* Constructors/Destructors                                                 */
-  /* ------------------------------------------------------------------------ */
-public:
-
-  ContactClusterCollection(ContactArea & area);
-  virtual ~ContactClusterCollection();
-
-  /* ------------------------------------------------------------------------ */
-  /* Methods                                                                  */
-  /* ------------------------------------------------------------------------ */
-public:
-
-  /// dump information about clusters to file
-  void dumpClusterInfoToFile(const std::string & filename);
-
-  Real getTotalArea();
-  Real getTotalPerimeter();
-  Real getNbClustersWithHole();
-  Real getNbClusters();
-
-
-  std::vector<Real> getAreas();
-  std::vector<Real> getMassCenters(UInt i);
-  std::vector<Real> getPerimeters();
-  std::vector<Real> getCompactnesses();
-  std::vector<Real> getInterPixelPerimeters();
-  std::vector<Real> getDiscretePerimeters();
-  std::vector<Real> getNbHoles();
-  //std::vector<Real> getRIs();
-  std::vector<Real> getColors();
-
-  /// write one file per cluster
-  void dumpAllClusters(const std::string & prefix);
-  /* ------------------------------------------------------------------------ */
-  /* Class Members                                                            */
-  /* ------------------------------------------------------------------------ */
-private:
-
-
-  //! the vector of clusters
-  std::vector<ContactCluster> clusters;
-
-  //! the input area where to extract all the clusters
-  ContactArea & area;
-
-  //! the total area
-  UInt TotalArea;
-  //! the total perimeter
-  UInt TotalPerimeter;
-  //! the number of clusters with holes
-  UInt NumberOfClustersWithHoles;
-
-};
-
-/* -------------------------------------------------------------------------- */
-
-__END_TAMAAS__
-
-#endif /* __CONTACT_CLUSTER_COLLECTION_HH__ */
diff --git a/src/grid.cpp b/src/grid.cpp
index ba1fc54..437f290 100644
--- a/src/grid.cpp
+++ b/src/grid.cpp
@@ -1,127 +1,128 @@
 /**
  *
  * @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 "tamaas.hh"
 #include "grid.hh"
 #include <cstring>
 #include <complex>
+#include <algorithm>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid():
   data() {
-  memset(this->n, 0, dim * sizeof(UInt));
-  memset(this->L, 0, dim * sizeof(Real));
+  std::fill(this->n, this->n+dim, 0);
+  std::fill(this->L, this->L+dim, 0);
   this->nb_components = 1;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim>::Grid(const UInt n[dim], const Real L[dim], UInt nb_components):
   data() {
-  memcpy(this->n, n, dim * sizeof(UInt));
-  memcpy(this->L, L, dim * sizeof(Real));
+  std::copy(n, n+dim, this->n);
+  std::copy(L, L+dim, this->L);
   this->nb_components = nb_components;
   this->resize(this->n);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 Grid<T, dim>::~Grid() {}
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::resize(const UInt *n) {
   if (n != this->n)
-    memcpy(this->n, n, dim * sizeof(UInt));
+    std::copy(n, n+dim, this->n);
   UInt size = this->computeSize();
   data.resize(size);
   data.assign(size, T(0.));
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::uniformSetComponents(const Grid<T, 1> &vec) {
   TAMAAS_ASSERT(vec.dataSize() == this->nb_components,
 		"Cannot set grid field with values of vector");
 
 #pragma omp parallel for
   for (UInt i = 0 ; i < dataSize() ; i++) {
     data[i] = vec(i % nb_components);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T, UInt dim>
 void Grid<T, dim>::printself(std::ostream & str) const {
   str << "Grid(" << dim << ", " << nb_components << ") {";
   for (UInt i = 0 ; i < data.size() - 1 ; i++) {
     str << data[i] << ", ";
   }
   str << data[data.size()-1] << "}";
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define GRID_SCALAR_OPERATOR_IMPL(op)			\
   template <typename T, UInt dim>			\
   inline void Grid<T, dim>::operator op(const T & e) {	\
     _Pragma("omp parallel for")				\
     for (UInt i = 0 ; i < this->data.size() ; i++) {	\
       data[i] op e;					\
     }							\
   }							\
 
 GRID_SCALAR_OPERATOR_IMPL(+=);
 GRID_SCALAR_OPERATOR_IMPL(*=);
 GRID_SCALAR_OPERATOR_IMPL(-=);
 GRID_SCALAR_OPERATOR_IMPL(/=);
 GRID_SCALAR_OPERATOR_IMPL(=);
 
 #undef GRID_SCALAR_OPERATOR_IMPL
 
 /* -------------------------------------------------------------------------- */
 
 /// Class instanciation
 #define GRID_INSTANCIATE_TYPE(type) template class Grid<type, 1>;	\
   template class Grid<type, 2>;						\
   template class Grid<type, 3>
 
 GRID_INSTANCIATE_TYPE(Real);
 GRID_INSTANCIATE_TYPE(UInt);
 GRID_INSTANCIATE_TYPE(Complex);
 GRID_INSTANCIATE_TYPE(int);
 GRID_INSTANCIATE_TYPE(bool);
 GRID_INSTANCIATE_TYPE(unsigned long);
 #undef GRID_INSTANCIATE_TYPE
 
 __END_TAMAAS__
diff --git a/src/isopowerlaw.cpp b/src/isopowerlaw.cpp
index 6ae78e3..1d9de73 100644
--- a/src/isopowerlaw.cpp
+++ b/src/isopowerlaw.cpp
@@ -1,57 +1,109 @@
 /**
  *
  * @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 "isopowerlaw.hh"
 #include "fftransform.hh"
+#include <map>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 template <UInt dim>
 void Isopowerlaw<dim>::computeFilter(GridHermitian<Real, dim> & filter_coefficients) const {
   const Real l[dim] = {1.};
   Grid<Real, dim> wavevectors(filter_coefficients.sizes(), l, dim);
   FFTransform<Real, dim>::template computeFrequencies<true>(wavevectors);
   auto end = wavevectors.end(dim);
 
 #pragma omp parallel for
   for (auto it = wavevectors.begin(dim) ; it < end ; ++it) {
     VectorProxy<Real> q(*it, dim);
     UInt i = wavevectors.getNbPoints() - (end - it);
     filter_coefficients(i) = (*this)(q);
   }
 }
 
+/* -------------------------------------------------------------------------- */
+
 template <UInt dim>
-Real Isopowerlaw<dim>::analyticalRMS() const {
+Real Isopowerlaw<dim>::rmsHeights() const {
   return std::sqrt(M_PI * ( (hurst + 1)/hurst * q1*q1
 			    - 1./hurst * std::pow(q1, 2*(hurst+1)) * std::pow(q2, -2*hurst)
 			    - q0*q0 ));
 }
 
+/* -------------------------------------------------------------------------- */
+
+/** 
+ * Analytical moments, cf. Yastrebov et al. (2015)
+ * "From infinitesimal to full contact between rough surfaces: Evolution
+ * of the contact area", appendix A
+ */
+template <>
+std::vector<Real> Isopowerlaw<2>::moments() const {
+  std::map<UInt, Real> T;
+  T[0] = 2*M_PI; T[2] = M_PI; T[4] = 3 * M_PI / 4.;
+  Real xi = q0 / q1;
+  Real zeta = q2 / q1;
+  
+  std::vector<Real> moments;
+
+  using namespace std;
+  for (int q : {0, 2, 4}) {
+    Real m = T[q] * pow(q1, q-2*hurst) * ((1-pow(xi, q+2))/(q+2) +
+					  (pow(zeta, q-2*hurst)-1)/(q-2*hurst));
+    moments.push_back(m);
+  }
+
+  return moments;
+}
+
+template <>
+std::vector<Real> Isopowerlaw<1>::moments() const {
+  TAMAAS_EXCEPTION("Moments have not been implemented for 1D surfaces");
+  return std::vector<Real>();
+}
+
+/* -------------------------------------------------------------------------- */
+
+template <UInt dim>
+Real Isopowerlaw<dim>::alpha() const {
+  std::vector<Real> m = moments();
+  return m[0] * m[2] / (m[1]*m[1]);
+}
+
+/* -------------------------------------------------------------------------- */
+
+template <UInt dim>
+Real Isopowerlaw<dim>::rmsSlopes() const {
+  std::vector<Real> m = moments();
+  return std::sqrt(M_PI * m[1] / 2);
+}
+
+
 template class Isopowerlaw<1>;
 template class Isopowerlaw<2>;
 
 __END_TAMAAS__
diff --git a/src/isopowerlaw.hh b/src/isopowerlaw.hh
index c5ad09e..412f1db 100644
--- a/src/isopowerlaw.hh
+++ b/src/isopowerlaw.hh
@@ -1,85 +1,98 @@
 /**
  *
  * @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 ISOPOWERLAW_H
 #define ISOPOWERLAW_H
 /* -------------------------------------------------------------------------- */
 #include "filter.hh"
 #include "grid_hermitian.hh"
 #include "types.hh"
 #include "tamaas.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /// Class representing an isotropic power law spectrum
 template <UInt dim>
 class Isopowerlaw : public Filter<dim> {
 public:
   /// Default constructor
   Isopowerlaw() {};
   /// Destructor
   virtual ~Isopowerlaw() {};
 
 public:
   /// Compute filter coefficients
   virtual void computeFilter(GridHermitian<Real, dim> & filter_coefficients) const;
 
   /// Compute a point of the PSD
   inline Real operator()(const VectorProxy<Real> & q_vec) const;
 
   /// Analytical rms of heights
-  Real analyticalRMS() const;
+  Real rmsHeights() const;
+
+  /// Analytical moments
+  std::vector<Real> moments() const;
+
+  /// Analytical Nayak's parameter
+  Real alpha() const;
+
+  /// Analytical RMS slopes
+  Real rmsSlopes() const;
 
 public:
   TAMAAS_ACCESSOR(q0, UInt, Q0);
   TAMAAS_ACCESSOR(q1, UInt, Q1);
   TAMAAS_ACCESSOR(q2, UInt, Q2);
   TAMAAS_ACCESSOR(hurst, Real, Hurst);
 
 protected:
   UInt q0, q1, q2;
   Real hurst;
 };
 
+/* -------------------------------------------------------------------------- */
+/* Inline implementations                                                     */
+/* -------------------------------------------------------------------------- */
+
 template<UInt dim>
 inline Real Isopowerlaw<dim>::operator()(const VectorProxy<Real> & q_vec) const {
   const Real C = 1.;
   const Real q = q_vec.l2norm();
 
   Real val;
 
   if (q < q0) val = 0;
   else if (q > q2) val = 0;
   else if (q < q1) val = C;
   else val = C * std::pow((q/q1), -2*(hurst+1));
 
   return std::sqrt(val);
 }
 
 __END_TAMAAS__
 
 #endif // ISOPOWERLAW_H
diff --git a/src/map_2d.cpp b/src/map_2d.cpp
deleted file mode 100644
index 1636479..0000000
--- a/src/map_2d.cpp
+++ /dev/null
@@ -1,141 +0,0 @@
-/**
- *
- * @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 "map_2d.hh"
-#include "fft_plan_manager.hh"
-#include <fstream>
-#include <sstream>
-#include <fftw3.h>
-#include <omp.h>
-/* -------------------------------------------------------------------------- */
-#define TIMER
-#include "surface_timer.hh"
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-
-template <typename T>
-void Map2d<T>::rescaleHeights(Real C){
-  /* compute actual RMS */
-  for (UInt i = 0 ; i < this->n[0] ; ++i)
-    for (UInt j = 0 ; j < this->n[1] ; ++j){
-      this->at(i,j) *= C;
-    }
-}
-
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-void Map2d<T>::setGridSize(UInt n0, UInt n1){
-  const UInt n[] = {n0, n1};
-  this->resize(n);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-void Map2d<T>::setGridSize(const UInt n[2]){
-  this->resize(n);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-void Map2d<T>::operator=(const T & e){
-  this->Grid<T, 2>::operator=(e);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-Map2d<T>::Map2d(UInt n0, UInt n1, Real L0, Real L1) {
-  setGridSize(n0,n1);
-  this->L[0] = L0;
-  this->L[1] = L1;
-}
-/* -------------------------------------------------------------------------- */
-template <typename T>
-Map2d<T>::Map2d(const UInt * n, const Real *L) {
-  setGridSize(n);
-  this->L[0] = L[0];
-  this->L[1] = L[1];
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-Map2d<T>::~Map2d()
-{}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-void Map2d<T>::dumpToTextFile(std::string filename) {
-
-  UInt n1 = this->size(0);
-  UInt n2 = this->size(1);
-  std::cout << "Writing surface file " << filename << std::endl;
-  std::ofstream file(filename.c_str());
-  if (!file.is_open()) SURFACE_FATAL("cannot open file " << filename);
-
-  file.precision(15);
-  for (UInt i = 0; i < n1; i++) {
-    for (UInt j=0; j< n2; j++) {
-      file << i << " " << j << " ";
-      this->writeValue(file,this->at(i,j));
-      file << std::endl;
-    }
-    file << std::endl;
-  }
-  file.close();
-}
-
-/* -------------------------------------------------------------------------- */
-template <typename T>
-void Map2d<T>::writeValue(std::ostream & os, const T & val) {
-  os << val;
-}
-
-/* -------------------------------------------------------------------------- */
-template <>
-void Map2d<Complex>::writeValue(std::ostream & os, const Complex & val) {
-
-  os << std::scientific << val.real();
-  os << " " << val.imag();
-}
-
-/* -------------------------------------------------------------------------- */
-
-template class Map2d<Complex>;
-template class Map2d<Real>;
-template class Map2d<int>;
-template class Map2d<unsigned int>;
-template class Map2d<bool>;
-template class Map2d<unsigned long>;
-
-__END_TAMAAS__
diff --git a/src/map_2d.hh b/src/map_2d.hh
deleted file mode 100644
index 41643ef..0000000
--- a/src/map_2d.hh
+++ /dev/null
@@ -1,416 +0,0 @@
-/**
- *
- * @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 __MAP_2D_HH__
-#define __MAP_2D_HH__
-/* -------------------------------------------------------------------------- */
-#include <vector>
-#include <iostream>
-#include <fftw3.h>
-#include "tamaas.hh"
-#include <fstream>
-#include "grid.hh"
-/* -------------------------------------------------------------------------- */
-
-__BEGIN_TAMAAS__
-
-template <typename T>
-class Map2d : public Grid<T, 2> {
-  /* ------------------------------------------------------------------------ */
-  /* Constructors/Destructors                                                 */
-  /* ------------------------------------------------------------------------ */
-public:
-  Map2d(UInt n1, UInt n2, Real L0 = 1., Real L1 = 1.);
-  Map2d(const UInt * n, const Real * L);
-
-  virtual ~Map2d();
-
-  /* ------------------------------------------------------------------------ */
-  /* Methods                                                                  */
-  /* ------------------------------------------------------------------------ */
-public:
-
-  template <bool periodic>
-  std::vector<std::pair<UInt,UInt> > getNeighborIndexes(UInt i, UInt j);
-  template <bool periodic>
-  std::vector<std::pair<UInt,UInt> > getNeighborIndexes(std::pair<UInt,UInt> i);
-
-  // //! return an iterator over the neighbors
-  // iterator beginNeighborhood();
-  // //! return the end iterator over the neighbors
-  // iterator endNeighborhood();
-
-  //compute compute component by component product and put result in object
-  template <typename T1, typename T2>
-  void fullProduct(Map2d<T1> & s1,Map2d<T2> &s2);
-  //! set the grid size
-  void setGridSize(const UInt n[2]);
-  //! set the grid size
-  void setGridSize(UInt n0,UInt n1);
-  //! rescale height profile
-  void rescaleHeights(Real factor);
-  //! Write heights array to text file usable with plot programs
-  void dumpToTextFile(std::string filename);
-  //! set to scalar value
-  void operator=(const T & val);
-  //! copy method
-  template <template<typename> class surf, typename T2>
-  void copy(const surf<T2> & s);
-  // copy the map passed as parameter
-  template <template<typename> class surf, typename T2>
-  void operator=(const surf<T2> & s);
-  // temporary hack
-  void operator=(const Map2d<Complex> & s);
-  // copy the map passed as const parameter
-  template <template<typename> class surf, typename T2>
-  void operator=(surf<T2> & s);
-
-/* -------------------------------------------------------------------------- */
-
-private:
-
-  //! write an entry to a stringstream
-  void writeValue(std::ostream & os, const T & val);
-
-  /* ------------------------------------------------------------------------ */
-  /* Accessors                                                                */
-  /* ------------------------------------------------------------------------ */
-
-public:
-
-  // access to the i,j-th value in the surface
-  inline T & operator()(UInt i,UInt j);
-  // access to the i,j-th value in the surface
-  inline T & operator()(std::pair<UInt,UInt> i);
-  // access to the i,j-th value in the surface
-  inline const T & operator()(UInt i,UInt j)const;
-  // access to the i,j-th value in the surface
-  inline const T & operator()(std::pair<UInt,UInt> i)const;
-  // access to the i,j-th value in the surface
-  inline T & at(UInt i,UInt j);
-  // access to the i-th value in the array
-  inline T & at(UInt i);
-  // access to the i,j-th value in the surface
-  inline T & at(std::pair<UInt,UInt> i);
-  // access to the i,j-th value in the surface
-  inline const T & at(UInt i,UInt j) const;
-  // access to the i,j-th value in the surface
-  inline const T & at(std::pair<UInt,UInt> i) const;
-  // access to the i,j-th value in the surface
-  inline const T & at(UInt i) const;
-  // access to the i-th value in the array
-  inline T & operator()(int i);
-  // access to the i-th value in the array
-  //  inline T & operator[](std::pair<UInt,UInt> i);
-  // access to the i-th value in the array
-  inline const T & operator()(int i) const;
-  //! get lateral number of discretization points
-  inline UInt size(UInt dir)const {return this->n[dir];};
-  //! get both lateral number of discretization points
-  inline const UInt* sizes()const {return this->n;};
-  //! get the lateral size of the sample
-  inline const Real * getLs() const {return this->L;};
-  //! get the lateral size of the sample
-  inline Real getL(UInt dir) const {return this->L[dir];};
-  //! get the lateral size of the sample
-  inline void setL(Real L[2]){this->L[0] = L[0];this->L[1] = L[1];};
-  //! get the lateral size of the sample
-  inline void setL(Real L0,Real L1){this->L[0] = L0;this->L[1] = L1;};
-
-  const T * getInternalData() const{return this->data.getPtr();};
-
-  //! Write heights array to text file usable with plot programs
-  template <bool consider_zeros=true>
-  void dumpToTextFile(std::string filename) const;
-};
-
-/* -------------------------------------------------------------------------- */
-//#define _at(I,J,N) I+J*N[1]
-#define _at(I,J,N) (I)*(N[1])+(J)
-/* -------------------------------------------------------------------------- */
-template <typename T>
-inline T & Map2d<T>::at(UInt i,UInt j){
-  return this->Grid<T, 2>::operator()(i, j, 0);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-inline const T & Map2d<T>::at(UInt i,UInt j)const {
-  return this->Grid<T, 2>::operator()(i, j, 0);
-}
-/* -------------------------------------------------------------------------- */
-#undef _at
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-inline T & Map2d<T>::at(UInt i){
-  return this->Grid<T, 2>::operator()(i);
-}
-
-/* -------------------------------------------------------------------------- */
-template <typename T>
-inline const T & Map2d<T>::at(UInt i) const{
-  return this->Grid<T, 2>::operator()(i);
-}
-
-/* -------------------------------------------------------------------------- */
-
-
-template <typename T>
-inline T & Map2d<T>::operator()(UInt i,UInt j){
-  return this->at(i,j);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-inline T & Map2d<T>::at(std::pair<UInt,UInt> i){
-  return this->at(i.first,i.second);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-inline const T & Map2d<T>::at(std::pair<UInt,UInt> i)const {
-  return this->at(i.first,i.second);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-inline T & Map2d<T>::operator()(std::pair<UInt,UInt> i){
-  return this->at(i.first,i.second);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-inline const T & Map2d<T>::operator()(UInt i,UInt j)const{
-  return this->at(i,j);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-inline const T & Map2d<T>::operator()(std::pair<UInt,UInt> i)const {
-  return this->at(i.first,i.second);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-inline T & Map2d<T>::operator()(int i){
-  return this->Grid<T, 2>::operator()(i);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-inline const T & Map2d<T>::operator()(int i) const{
-  return this->Grid<T, 2>::operator()(i);
-}
-
-/* -------------------------------------------------------------------------- */
-// template <typename T>
-// inline T & Map2d<T>::operator[](std::pair<UInt,UInt> c){
-//   return this->at(c.first,c.second);
-// }
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-template <bool periodic>
-std::vector<std::pair<UInt,UInt> > Map2d<T>::getNeighborIndexes(UInt i, UInt j){
-
-  std::vector<std::pair<UInt,UInt> > vres;
-
-
-  std::pair<int,int> res;
-
-  for (int ii = -1; ii < 2; ++ii) {
-    res.first = ii+static_cast<int>(i);
-    if (res.first < 0 && !periodic)    continue;
-    if (res.first >= (int)this->n[0] && !periodic) continue;
-    if (res.first < 0 && periodic)     res.first += this->n[0];
-    if (res.first >= (int)this->n[0] && periodic) res.first -= this->n[0];
-
-    for (int jj = -1; jj < 2; ++jj) {
-
-    res.second = jj+static_cast<int>(j);
-    if (res.second < 0 && !periodic)    continue;
-    if (res.second >= (int)this->n[1] && !periodic) continue;
-    if (res.second < 0 && periodic)     res.second += this->n[1];
-    if (res.second >= (int)this->n[1] && periodic) res.second -= this->n[1];
-
-    if (ii == 0 && jj == 0) continue;
-
-    {
-      UInt i = res.first;
-      UInt j = res.second;
-
-      if (i+j*this->n[0] >= this->n[0]*this->n[1]) SURFACE_FATAL("out of bound access "
-					     << i << " " << j << " "
-					     << this->n[0] << " " << this->n[1] << " "
-					     << i+j*this->n[0]);
-    }
-    vres.push_back(res);
-    }
-  }
-  return vres;
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-template <bool periodic>
-std::vector<std::pair<UInt,UInt> > Map2d<T>
-::getNeighborIndexes(std::pair<UInt,UInt> i){
-  return this->getNeighborIndexes<periodic>(i.first,i.second);
-}
-
-/* -------------------------------------------------------------------------- */
-template <typename T>
-template <typename T1, typename T2>
-void Map2d<T>::fullProduct(Map2d<T1> & s1,Map2d<T2> & s2){
-
-  auto n1 = s1.sizes();
-  auto n2 = s2.sizes();
-
-  if (n1[0] != n2[0] ||
-      n1[1] != n2[1]) SURFACE_FATAL("fullproduct of non compatible matrices");
-
-  if (this->n[0] != n1[0] ||
-      this->n[1] != n1[1]) SURFACE_FATAL("fullproduct of non compatible matrices")
-
-  for (UInt i = 0 ; i < this->n[0] ; ++i)
-    for (UInt j = 0 ; j < this->n[1] ; ++j){
-      this->at(i,j) = s1(i,j)* s2(i,j);
-    }
-}
-/* -------------------------------------------------------------------------- */
-template <typename T>
-template <bool consider_zeros>
-void Map2d<T>::dumpToTextFile(std::string filename) const {
-
-  UInt m = this->n[0];
-  UInt n = this->n[1];
-  std::cout << "Writing surface file " << filename << std::endl;
-  std::ofstream file(filename.c_str());
-  if (!file.is_open()) SURFACE_FATAL("cannot open file " << filename);
-
-  file.precision(15);
-  for (UInt i = 0; i < m; i++) {
-    for (UInt j=0; j< n; j++) {
-      if (!consider_zeros &&
-	  std::norm(this->at(i,j)) < zero_threshold*zero_threshold)
-	continue;
-
-      file << i << " " << j << " ";
-      file << std::scientific << this->at(i,j);
-      file << std::endl;
-    }
-    file << std::endl;
-  }
-  file.close();
-}
-
-
-// template <typename T>
-// class Map2d<T>::iterator {
-
-// public:
-
-//   iterator(std::vector<T> & data,
-// 	   std::queue<std::pair<UInt,UInt> > & indices):
-//     data(data),indices(indices){
-//   }
-
-//   inline void operator++(){
-//     indices.pop();
-//   }
-
-//   inline T & operator *(){
-//     return data[indices.front];
-//   }
-
-//   bool operator != (Map2d<T>::iterator & it){
-//     return (indices.front == it.indices.front);
-//   }
-
-// private:
-
-//   std::vector<T> & data;
-//   std::queue<std::pair<UInt,UInt> >  indices;
-
-//   bool finished;
-// };
-
-
-/* -------------------------------------------------------------------------- */
-template <typename T>
-template <template<typename> class surf, typename T2>
-void Map2d<T>::copy(const surf<T2> & s){
-  //this->setGridSize(s.sizes());
-  //UInt m = this->sizes()[0];
-  //UInt n = this->sizes()[1];
-  //#pragma omp parallel for
-  //for (UInt i = 0; i < m*n; i++) {
-  //  this->at(i) = s(i);
-  //}
-  this->Grid<T, 2>::copy(s);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-template <template<typename> class surf, typename T2>
-void Map2d<T>::operator=(surf<T2> & s){
-  this->copy(s);
-}
-/* -------------------------------------------------------------------------- */
-template <typename T>
-template <template<typename> class surf, typename T2>
-void Map2d<T>::operator=(const surf<T2> & s){
-  this->copy(s);
-}
-
-/* -------------------------------------------------------------------------- */
-
-template <typename T>
-void Map2d<T>::operator=(const Map2d<Complex> &){
-  SURFACE_FATAL("Cannot copy a complex surface into incompatible type");
-}
-template <>
-inline void Map2d<Complex>::operator=(const Map2d<Complex> & s){
-  this->copy(s);
-}
-
-/* -------------------------------------------------------------------------- */
-
-__END_TAMAAS__
-
-#endif /* __MAP_2D_HH__ */
diff --git a/src/map_2d_square.cpp b/src/map_2d_square.cpp
index 0c9dac7..e1eaece 100644
--- a/src/map_2d_square.cpp
+++ b/src/map_2d_square.cpp
@@ -1,106 +1,106 @@
 /**
  *
  * @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 "map_2d_square.hh"
 #include <fstream>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 
 template <typename T>
 void Map2dSquare<T>::loadFromTextFile(std::string filename, int each) {
 
   std::string line;
   std::ifstream myfile (filename.c_str());
   UInt n = 0;
   if (myfile.is_open()) {
       while (! myfile.eof() ) {
 
 	  getline (myfile,line);
 	  if (line[0] == '#')
 	    continue;
 	  if (line.size() == 0)
 	    continue;
 	  std::stringstream s;
 	  s << line;
 	  UInt i,j;
 	  s >> i;
 	  s >> j;
 	  if (n < i) n = i;
 	  if (n < j) n = j;
 	}
       ++n;
       std::cout << "read " << n << " x " << n << " grid points" << std::endl;
 
       myfile.clear();              // forget we hit the end of file
       myfile.seekg(0, std::ios::beg);   // move to the start of the file
 
       n = int(n/Real(each));
       setGridSize(n);
 
       while (! myfile.eof() ) {
 	getline (myfile,line);
 	if (line == "") continue;
 	std::stringstream s;
 	s << line;
 	int i,j;
 	s >> i;
 	s >> j;
 	T val = this->readValue(s);
 	if (i % each == 0 && j % each == 0) {
-	  this->at(i/each,j/each) = val;
+	  (*this)(i/each,j/each) = val;
 	}
       }
       myfile.close();
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 T Map2dSquare<T>::readValue(std::stringstream & sstr) {
   SURFACE_FATAL("To implement");
 }
 
 /* -------------------------------------------------------------------------- */
 template <>
 Complex Map2dSquare<Complex>::readValue(std::stringstream & sstr) {
 
   Real h = 0.,hIm = 0.;
   sstr >> h;
   sstr >> hIm;
   return Complex(h,hIm);
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 template class Map2dSquare<Complex>;
 template class Map2dSquare<Real>;
 template class Map2dSquare<int>;
 template class Map2dSquare<unsigned long>;
 
 __END_TAMAAS__
diff --git a/src/map_2d_square.hh b/src/map_2d_square.hh
index 37282c1..3f6eb22 100644
--- a/src/map_2d_square.hh
+++ b/src/map_2d_square.hh
@@ -1,81 +1,93 @@
 /**
  *
  * @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 __MAP_2D_SQUARE_HH__
 #define __MAP_2D_SQUARE_HH__
 /* -------------------------------------------------------------------------- */
-#include "map_2d.hh"
+#include "grid.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 template <typename T>
-class Map2dSquare : public Map2d<T> {
+class Map2dSquare : public Grid<T, 2> {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 public:
 
-  Map2dSquare(UInt n, Real L = 1.): Map2d<T>(n,n,L,L){};
-  virtual ~Map2dSquare(){};
+  Map2dSquare(UInt n, Real L = 1.): Grid<T, 2>() {
+    UInt dims[2] = {n, n};
+    this->resize(dims);
+    this->L[0] = L;
+    this->L[1] = L;
+  }
+
+  virtual ~Map2dSquare() {}
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 public:
 
   //! get lateral number of discretization points
-  inline UInt size() const{return Map2d<T>::size(0);};
-  //! load all n-th point of heights array from text file usable with plot programs
-  void loadFromTextFile(std::string filename, int each=1);
+  inline UInt size() const{return this->n[0];};
 
   //! set the grid size
-  void setGridSize(UInt n){Map2d<T>::setGridSize(n,n);};
+  void setGridSize(UInt n) {
+    UInt dims[2] = {n, n};
+    this->resize(dims);
+  }
+
+  void loadFromTextFile(std::string filename, int each);
 
 
 private:
   //! read an entry from a stringstream
   T readValue(std::stringstream & sstr);
 
   /* ------------------------------------------------------------------------ */
   /* Accessors                                                                */
   /* ------------------------------------------------------------------------ */
 public:
 
   //! get the lateral size of the sample
-  inline Real getL()const {return Map2d<T>::getL(0);};
+  inline Real getL() const {return this->L[0];};
   //! get the lateral size of the sample
-  inline void setL(Real L){Map2d<T>::setL(L,L);};
+  inline void setL(Real L) {
+    this->L[0] = L;
+    this->L[1] = L;
+  }
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
 
 };
 
 __END_TAMAAS__
 #endif /* __MAP_2D_SQUARE_HH__ */
diff --git a/src/surface.cpp b/src/surface.cpp
index 99c3547..7aaca73 100644
--- a/src/surface.cpp
+++ b/src/surface.cpp
@@ -1,445 +1,445 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 #include <iostream>
 #include <fstream>
 #include <sstream>
 #include <string>
 #include <cmath>
 #include <random>
 /* -------------------------------------------------------------------------- */
 #include "surface.hh"
 #include "surface_statistics.hh"
 /* -------------------------------------------------------------------------- */
 #ifdef USING_IOHELPER
 #include "dumper_paraview.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 std::string fname = "surface_height.pdf";
 /* -------------------------------------------------------------------------- */
 
 
 
 template <typename T>
 void Surface<T>::expandByLinearInterpolation() {
   UInt n = this->size();
   Surface * old_s = new Surface(n,1.);
   int new_n = n*2;
   *old_s = *this;
   this->setGridSize(new_n);
   for (UInt i=0; i<n; i++) {
     for (UInt j=0; j<n; j++) {
       (*this)(i,j) = (*old_s)(i,j);
     }
   }
   delete old_s;
 }
 
 Real sinc(Real x){
   return sin(M_PI*x)/(M_PI*x);
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByShannonInterpolation(int factor) {
 
   /* /!\ this->FFTTransform(); /!\ */
   TAMAAS_EXCEPTION("Shannon interpolation needs old FFTTransform interface");
   UInt n = this->size();
   UInt N = std::pow(2,factor)*this->size();
   std::cout << N << " " << factor << std::endl;
   Surface s = Surface(N,1.);
 
 
   for (UInt i = 0; i < n/2; ++i)
     for (UInt j = 0; j < n/2; ++j) {
-      s(i,j) = this->at(i,j);
+      s(i,j) = (*this)(i,j);
     }
 
   for (UInt i = n/2; i < n; ++i)
     for (UInt j = 0; j < n/2; ++j) {
-      s(N-n+i,j) = this->at(i,j);
+      s(N-n+i,j) = (*this)(i,j);
     }
 
   for (UInt i = n/2; i < n; ++i)
     for (UInt j = n/2; j < n; ++j) {
-      s(N-n+i,N-n+j) = this->at(i,j);
+      s(N-n+i,N-n+j) = (*this)(i,j);
     }
 
   for (UInt i = 0; i < n/2; ++i)
     for (UInt j = n/2; j < n; ++j) {
-      s(i,N-n+j) = this->at(i,j);
+      s(i,N-n+j) = (*this)(i,j);
     }
 
   //  s.dumpToTextFile<false,true>("temp.csv");
 
   /* /!\ s.FFTITransform(); /!\ */
 
   this->setGridSize(N);
   Real f = std::pow(2,factor);
   f = f*f;
 
   for (UInt i = 0; i < N*N; ++i)
-    this->at(i) = f*s(i);
+    (*this)(i) = f*s(i);
 
 
 //   int new_n = 2*n;
 //   Surface s(*this);
 //   this->setGridSize(new_n);
 
 //   // put the original signal with zeros for to be interpolated points
 //   for (int i=0; i<n; i++)
 //     for (int j=0; j<n; j++)
-//       this->at(2*i,j*2) = s(i,j);
+//       (*this)(2*i,j*2) = s(i,j);
 
 //   // do shannon interpolation
 //   int rcut = 64;
 
 //   //interpolate lines
 //   for (int i=0; i<new_n; i++) {
 //     for (int j=0; j<new_n; j++) {
 //       if (i%2 == 0) continue;
 //       // Real loop on local neighbors to interpolate
 //       for (int k=-rcut; k<rcut+1; k++) {
 // 	int i2 = i + k;
 // 	if (i2 < 0) i2 += new_n;
 // 	if (i2 >= new_n) i2 -= new_n;
 // 	if (i2%2) continue;
 
 // 	Real d = fabs(k);
 // 	s[i+j*new_n] += s(i2,j)*sinc(d/2);
 //       }
 //     }
 //   }
 
 //   //interpolate columns
 //   for (int i=0; i<new_n; i++) {
 //     for (int j=0; j<new_n; j++) {
 //       if (j%2 == 0) continue;
 //       // Real loop on local neighbors to interpolate
 //       for (int k=-rcut; k<rcut+1; k++) {
 // 	int j2 = j + k;
 // 	if (j2 < 0) j2 += new_n;
 // 	if (j2 >= new_n) j2 -= new_n;
 // 	if (j2%2) continue;
 
 // 	Real d = fabs(k);
 // 	s[i+j*new_n] += complex(s[i+j2*new_n].real()*sinc(d/2),
 // 				s[i+j2*new_n].imag()*sinc(d/2));
 //       }
 //     }
 //   }
 
 
 // //   for (int i=0; i<n; i++) {
 // //     for (int j=0; j<n; j++) {
 // //       s[2*i+j*2*new_n].real() = 0;
 // //       s[2*i+j*2*new_n].im = 0;
 // //     }
 // //   }
 
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByMirroring() {
   UInt n = this->size();
   UInt new_n = n*2;
   Surface s(*this);
   this->setGridSize(new_n);
 
   for (UInt i=0; i<new_n; i++) {
     for (UInt j=0; j<new_n; j++) {
       UInt i1 = i;
       if (i1 > n) i1 = new_n - i1;
       UInt j1 = j;
       if (j1 > n) j1 = new_n - j1;
       if (i == n || j == n)  {
-	this->at(i,j) = 0.;
+	(*this)(i,j) = 0.;
 	continue;
       }
-      this->at(i,j) = s(i1,j1);
+      (*this)(i,j) = s(i1,j1);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::expandByPuttingZeros() {
   UInt n = this->size();
   UInt new_n = n*2;
   Surface s(*this);
   this->setGridSize(new_n);
 
   for (UInt i=0; i< n; i++) {
     for (UInt j=0; j< n; j++) {
-      this->at(2*i,j*2) = s(i,j);
+      (*this)(2*i,j*2) = s(i,j);
     }
   }
 }
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::contractByTakingSubRegion() {
   UInt n = this->size();
   Surface s(*this);
   this->setGridSize(n/2);
 
   for (UInt i=0; i < n/2; i++) {
     for (UInt j=0; j < n/2; j++) {
-      this->at(i,j) = s(i,j);
+      (*this)(i,j) = s(i,j);
     }
   }
 }
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::contractByIgnoringMidPoints() {
   UInt n = this->size();
   Surface s(*this);
   this->setGridSize(n/2);
 
   for (UInt i=0; i<n; i++) {
     for (UInt j=0; j<n; j++) {
       if (i%2) continue;
       if (j%2) continue;
-      this->at(i/2,j/2) = s(i,j);
+      (*this)(i/2,j/2) = s(i,j);
     }
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 
 // void Surface<T>::dumpToTextFile(string filename, bool pr) {
 
 //   cout << "Writing surface file " << filename << endl;
 //   ofstream file(filename.c_str());
 //   file.precision(15);
 
 //   if (pr)
 //   {
 //     for (int i=0; i<n; i++) {
 //       for (int j=0; j<n; j++) {
 //         file << i << " " << j << " ";
 //         if (heights[i+j*n].real() < 1e-10)
 //           file << scientific << -1e10 << " " << 0;
 //         else
 //           file << scientific << heights[i+j*n].real() << " " << heights[i+j*n].im;
 //         file << endl;
 //       }
 //       file << endl;
 //     }
 //   }
 //   else
 //   {
 //     for (int i=0; i<n; i++) {
 //       for (int j=0; j<n; j++) {
 //         file << i << " " << j << " ";
 //         file << scientific << heights[i+j*n].real() << " " << heights[i+j*n].im;
 //         file << endl;
 //       }
 //     }
 //   }
 //   file.close();
 
 // }
 
 /* -------------------------------------------------------------------------- */
 
 
 
 // Surface * Surface<T>::loadFromTextFile(std::string filename) {
 
   // std::string line;
   // std::ifstream myfile (filename.c_str());
 
   // if (!myfile.is_open()){
   //   SURFACE_FATAL("cannot open file " << filename);
   // }
 
   // Surface * my_surface = NULL;
   // UInt n = 0;
   // if (myfile.is_open())
   //   {
   //     while (! myfile.eof() )
   // 	{
   // 	  getline (myfile,line);
   // 	  if (line[0] == '#')
   // 	    continue;
   // 	  if (line.size() == 0)
   // 	    continue;
 
   // 	  std::stringstream s;
   // 	  s << line;
   // 	  int i,j;
   // 	  s >> i;
   // 	  s >> j;
   // 	  if (n < i) n = i;
   // 	  if (n < j) n = j;
   // 	}
   //     ++n;
   //     std::cout << "read " << n << " x " << n << " grid points" << std::endl;
 
   //     myfile.clear();              // forget we hit the end of file
   //     myfile.seekg(0, std::ios::beg);   // move to the start of the file
 
   //     my_surface = new Surface(n);
 
   //     while (! myfile.eof() )
   // 	{
   // 	  getline (myfile,line);
   // 	  if (line[0] == '#')
   // 	    continue;
   // 	  if (line.size() == 0)
   // 	    continue;
 
   // 	  std::stringstream s;
   // 	  s << line;
   // 	  int i,j;
   // 	  Real h,hIm;
   // 	  s >> i;
   // 	  s >> j;
   // 	  s >> h;
   // 	  s >> hIm;
   // 	  my_surface->heights[i+j*n] = complex(h,hIm);
   // 	}
   //     myfile.close();
   //   }
 //   return my_surface;
 // }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void Surface<T>::dumpToParaview(std::string filename) const {
 #ifdef USING_IOHELPER
 
   UInt n = this->size();
   const Surface & data = *this;
   /* build 3D coordinates array to be usable with iohelper */
   //  cerr << "allocate coords for paraview" << endl;
   Real *  coords = new Real[n*n*3];
   //Real *  coords = new Real[(2*n-1)*(2*n-1)*3];
   //cerr << "allocate zcoords for paraview" << endl;
   Real *  zcoords = new Real[n*n*3];
   //Real *  zcoords = new Real[(2*n-1)*(2*n-1)*3];
   //cerr << "allocate zcoords2 for paraview" << endl;
   Real *  zcoords2 = new Real[n*n*3];
   //  Real *  zcoords2 = new Real[(2*n-1)*(2*n-1)*3];
   int index = 0;
   // for (int i=-n+1; i<n; i++) {
   //   for (int j=-n+1; j<n; j++) {
   for (int i=0; i<n; i++) {
     for (int j=0; j<n; j++) {
       coords[3*index] = 1.*i/n;
       coords[3*index+1] = 1.*j/n;
       coords[3*index+2] = 0;
       zcoords[3*index] = 0;
       zcoords[3*index+1] = 0;
       zcoords2[3*index] = 0;
       zcoords2[3*index+1] = 0;
       int i1 = (i%n);
       int j1 = (j%n);
       if (i1 < 0) i1 += n;
       if (j1 < 0) j1 += n;
       zcoords[3*index+2] = data(i1,j1).real();
       zcoords2[3*index+2] = data(i1,j1).imag();
       ++index;
     }
   }
 
   iohelper::DumperParaview dumper;
   dumper.setPoints(coords, 3, index, filename.c_str());
   dumper.setConnectivity(NULL, iohelper::POINT_SET, 1, iohelper::FORTRAN_MODE);
   dumper.addNodeDataField("data.re",zcoords,3,index);
   dumper.addNodeDataField("data.im",zcoords2,3,index);
   dumper.setPrefix("./paraview");
   dumper.init();
   dumper.dump();
 
   delete [] coords;
   delete [] zcoords;
   delete [] zcoords2;
 #endif //USING_IOHELPER
 }
 /* -------------------------------------------------------------------------- */
 
 
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 Surface<T>::Surface(UInt a, Real L) : Map2dSquare<T>(a,L){
   this->setGridSize(a);
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 Surface<T>::~Surface(){
 }
 
 /* -------------------------------------------------------------------------- */
 
 // void Surface<T>::copy(Surface & s){
 //   this->n = s.n;
 //   this->L = s.L;
 //   this->heights = s.heights;
 // }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::recenterHeights(){
   UInt n = this->size();
   Real zmean = SurfaceStatistics::computeAverage(*this);
 
   for (UInt i = 0 ; i < n*n ; ++i)
-    this->at(i) = this->at(i) - zmean;
+    (*this)(i) = (*this)(i) - zmean;
 }
 
 /* -------------------------------------------------------------------------- */
 template <typename T>
 void Surface<T>::generateWhiteNoiseSurface(Surface<Real> & sources, long random_seed){
   //set the random gaussian generator
   std::mt19937 gen(random_seed);
   std::normal_distribution<> gaussian_generator;
 
   UInt n = sources.size();
 
   for (UInt i = 0 ; i < n ; ++i)
     for (UInt j = 0 ; j < n ; ++j){
       sources(i,j) = gaussian_generator(gen);
     }
 }
 /* -------------------------------------------------------------------------- */
 
 template class Surface<Real>;
 template class Surface<unsigned long>;
 template class Surface<UInt>;
 template class Surface<int>;
 template class Surface<bool>;
 
 __END_TAMAAS__
diff --git a/src/surface.hh b/src/surface.hh
index 5ef47fa..d2f49b3 100644
--- a/src/surface.hh
+++ b/src/surface.hh
@@ -1,92 +1,92 @@
 /**
  *
  * @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_HH__
 #define __SURFACE_HH__
 /* -------------------------------------------------------------------------- */
 #include <string>
 #include <math.h>
 #include <iostream>
 #include <limits>
 #include <complex>
 /* -------------------------------------------------------------------------- */
 #include "map_2d_square.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 template <typename T>
 class Surface : public Map2dSquare<T> {
 
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 
  public:
 
   Surface(UInt a, Real L);
   ~Surface();
 
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 
   //! Write VTU file that is a set of points
   void dumpToParaview(std::string filename) const;
   //! create a new surface with linear interpolated points
   void expandByLinearInterpolation();
   //! create a new surface with linear interpolated points
   void expandByShannonInterpolation(int factor);
   //! create a new surface with mirroring
   void expandByMirroring();
   //! create a new surface with mirroring
   void expandByPuttingZeros();
   //! contract the surface by ignoring extra data
   void contractByIgnoringMidPoints();
   //! contract the surface by taking left bottom part of it
   void contractByTakingSubRegion();
   //recenter the heights to zero value
   void recenterHeights();
   // set all the members of the map to a given value
-  using Map2d<T>::operator=;
+  using Grid<T, 2>::operator=;
 
 
   static void generateWhiteNoiseSurface(Surface<Real> & surface, long random_seed);
 
   virtual const Surface<T> & getWrappedNumpy(){return *this;};
 
 };
 
 /* -------------------------------------------------------------------------- */
 
 /* -------------------------------------------------------------------------- */
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 
 #include "surface_complex.hh"
 
 #endif /* __SURFACE_HH__ */
diff --git a/src/surface_generator.hh b/src/surface_generator.hh
index 09329ea..07332bd 100644
--- a/src/surface_generator.hh
+++ b/src/surface_generator.hh
@@ -1,64 +1,64 @@
 /**
  *
  * @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_GENERATOR_H
 #define SURFACE_GENERATOR_H
 /* -------------------------------------------------------------------------- */
 #include "grid.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /// Class generating random surfaces
 template<UInt dim>
 class SurfaceGenerator {
 public:
   /// Default constructor
   SurfaceGenerator();
   /// Default destructor
   virtual ~SurfaceGenerator();
 
 public:
   /// Build surface profile (array of heights)
   virtual Grid<Real, dim> & buildSurface() = 0;
   /// Set surface sizes
   void setSizes(const UInt n[dim]);
-  /// Get random seed
-  long & getRandomSeed() { return random_seed; }
+
+  TAMAAS_ACCESSOR(random_seed, long, RandomSeed);
 
   /// Maintaining old interface
   friend class SurfaceGeneratorFilterFFT;
 
 protected:
   Grid<Real, dim> grid;
   long random_seed;
 };
 
 /* -------------------------------------------------------------------------- */
 
 __END_TAMAAS__
 
 #endif
diff --git a/src/surface_generator_filter_fft.hh b/src/surface_generator_filter_fft.hh
index 099d565..3febbb8 100644
--- a/src/surface_generator_filter_fft.hh
+++ b/src/surface_generator_filter_fft.hh
@@ -1,80 +1,80 @@
 /**
  *
  * @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_GENERATOR_FILTER_FFT_H
 #define SURFACE_GENERATOR_FILTER_FFT_H
 /* -------------------------------------------------------------------------- */
 #include "surface_generator_filter.hh"
 #include "isopowerlaw.hh"
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 typedef Grid<Real, 2> Grid2dReal;
 
 /// Legacy class to maintain interface
 class SurfaceGeneratorFilterFFT {
  public:
   //! can construct a generator with default values
   SurfaceGeneratorFilterFFT();
   //! destructor
   virtual ~SurfaceGeneratorFilterFFT();
 
 public:
   //! get high frequency cut
   UInt & getQ2(){return pl_filter.getQ2();}
   //! get low frequency cut
   UInt & getQ0(){return pl_filter.getQ0();}
   //! get roll frequency cut
   UInt & getQ1(){return pl_filter.getQ1();}
   //! get hurst
   Real & getHurst(){return pl_filter.getHurst();}
   /// compute analytical RMS
-  Real analyticalRMS() const {return pl_filter.analyticalRMS(); }
+  Real analyticalRMS() const {return pl_filter.rmsHeights(); }
   /// Swig does not wrap the numpy otherwise
   Grid2dReal & buildSurface();
   /// Get gridsize
   UInt & getGridSize() { return grid_size; }
   /// Get random seed
   long & getRandomSeed() { return generator.getRandomSeed(); }
   /// Get RMS
   Real & getRMS() { return rms; }
   /// Contrain RMS
   Real constrainRMS();
   /// Init
   void Init() {}
 
 protected:
   SurfaceGeneratorFilter<2> generator;
   Isopowerlaw<2> pl_filter;
 
   UInt grid_size;
   Real rms;
 };
 
 __END_TAMAAS__
 
 #endif
diff --git a/src/tamaas.hh b/src/tamaas.hh
index 643a99c..b254164 100644
--- a/src/tamaas.hh
+++ b/src/tamaas.hh
@@ -1,105 +1,105 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  * @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 TAMAAS_HH
 #define TAMAAS_HH
 /* -------------------------------------------------------------------------- */
 #include <exception>
 #include <string>
 #include <complex>
 #include <iostream>
 /* -------------------------------------------------------------------------- */
 
 /// Namespace macros
 #define __BEGIN_TAMAAS__ namespace tamaas {
 #define __END_TAMAAS__ }
 
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 /// Common types definitions
 typedef double Real;
 typedef unsigned int UInt;
 typedef int Int;
 typedef std::complex<Real> Complex;
 static const Real zero_threshold = 1e-14;
 
 /* -------------------------------------------------------------------------- */
 /// initialize tamaas (0 threads => let OMP_NUM_THREADS decide)
 void initialize(UInt num_threads = 0);
 /// cleanup tamaas
 void finalize();
 /* -------------------------------------------------------------------------- */
 class Exception : public std::exception {
 public:
   Exception(const std::string & mesg) throw() {
     msg = mesg;
   }
 
   virtual const char* what() const throw() {
     return msg.c_str();
   }
 
   virtual ~Exception() throw() {}
 
 private:
   std::string msg;
 };
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 #define TAMAAS_EXCEPTION(mesg)	{			\
     std::stringstream sstr;				\
     sstr						\
       <<  __FILE__					\
       << ":" << __LINE__ << ":FATAL: "			\
       <<  mesg << std::endl;				\
     std::cerr.flush();					\
     throw ::tamaas::Exception(sstr.str());		\
   }
 
 #define SURFACE_FATAL(mesg) TAMAAS_EXCEPTION(mesg)
 
 #if defined(TAMAAS_DEBUG)
 #define TAMAAS_ASSERT(cond, reason) do {			\
     if (!(cond)) {						\
       TAMAAS_EXCEPTION(#cond " assert failed: " << reason);	\
     } } while(0)
 #define TAMAAS_DEBUG_EXCEPTION(reason) TAMAAS_EXCEPTION(reason)
 #else
-#define TAMAAS_ASSERT(mesg, reason)
+#define TAMAAS_ASSERT(cond, reason)
 #define TAMAAS_DEBUG_EXCEPTION(reason)
 #endif
 #define TAMAAS_ACCESSOR(var, type, name) type & get##name() { return var; }	\
   void set##name(const type & new_var) { var = new_var; }
 
 /* -------------------------------------------------------------------------- */
 
 #endif // TAMAAS_HH
diff --git a/src/types.hh b/src/types.hh
index bba9353..f8ac28a 100644
--- a/src/types.hh
+++ b/src/types.hh
@@ -1,106 +1,106 @@
 /**
  *
  * @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 __TYPES_HH__
 #define __TYPES_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid.hh"
 #include <cstddef>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 template <typename T, UInt dim>
 class Proxy : public Grid<T, dim> {
 public:
   /// Default constructor
   Proxy() : Grid<T, dim>() {}
   /// Copy constructor
   Proxy(Proxy & p) : Grid<T, dim>() {
     if (p.n != this->n)
-      memcpy(this->n, p.n, dim * sizeof(UInt));
+      std::copy(p.n, p.n+dim, this->n);
     this->nb_components = p.nb_components;
     this->data.wrapMemory(p.getInternalData(), p.dataSize());
   }
   /// Destructor
   virtual ~Proxy() {}
 public:
   /// Set pointer
   void setPointer(T * ptr) {this->data.setPtr(ptr);}
   /// Get increment
   ptrdiff_t increment() const {return this->data.size();}
 };
 
 // Forward declarations
 template <typename T> class MatrixProxy;
 
 template <typename T>
 class VectorProxy : public Proxy<T, 1> {
 public:
   /// Wrap on memory
   VectorProxy(T * mem, UInt n);
   /// Copy constructor
   VectorProxy(VectorProxy & v): Proxy<T, 1>(v) {}
   /// Destructor
   virtual ~VectorProxy();
 
 public:
   /// Access
   const T & operator()(UInt i) const;
   /// Access (non-const)
   T & operator()(UInt i);
   /// Matrix-vector multiplication
   void mul(const MatrixProxy<T> & mat, const VectorProxy<T> & vec);
   /// Dot product
   T dot(const VectorProxy<T> & vec);
   /// L2 norm
   T l2norm() const;
   /// = operator
   using Grid<T, 1>::operator=;
 };
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 class MatrixProxy : public Proxy<T, 2> {
 public:
   /// Wrap on memory
   MatrixProxy(T * mem, UInt n, UInt m);
   /// Copy constructor
   MatrixProxy(MatrixProxy & v): Proxy<T, 2>(v) {}
   /// Destructor
   virtual ~MatrixProxy();
 
 public:
   /// Access
   const T & operator()(UInt i, UInt j) const;
   /// Access (non-const)
   T & operator()(UInt i, UInt j);
 };
 
 __END_TAMAAS__
 
 #endif // __TYPES_HH__
diff --git a/tests/test_surface.py b/tests/test_surface.py
index f7efd80..8fe12ca 100644
--- a/tests/test_surface.py
+++ b/tests/test_surface.py
@@ -1,147 +1,135 @@
 #!/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/>.
 # -----------------------------------------------------------------------------
 from __future__ import print_function, division
 
 import sys
 import numpy as np
 import tamaas as tm
 
 
 def compute_mse(values, mean):
     return np.mean((values - mean)**2) / values.size
 
 
 def main():
     tm.initialize(1)
 
     # Spectrum parameters in wavelengths
     disctretization_size = 1.
     lambda_s = 8.
     lambda_r = 64.
     lambda_l = 64.
     box_size = 512
 
     # Independent parameters
     Hurst = 0.8
     hrms = 1.
 
     # Spectrum parameters in wavenumbers
     k_l = int(box_size // lambda_l)
     k_r = int(box_size // lambda_r)
     k_s = int(box_size // lambda_s)
     N = int(box_size // disctretization_size)
 
-    # Analytical moments, cf. Yastrebov et al. (2015)
-    # "From infinitesimal to full contact between rough surfaces: Evolution
-    # of the contact area", appendix A
-    xi = k_l / k_r
-    zeta = k_s / k_r
-    H = Hurst
-    T = {0: 2*np.pi,
-         2: np.pi,
-         4: 3*np.pi/4}
-
-    m = {}
-    for q in [0, 2, 4]:
-        m[q] = T[q] * k_r**(q-2*H)*((1-xi**(q+2))/(q+2) + (zeta**(q-2*H)-1)/(q-2*H))
-
-    # Some spectral parameters
-    alpha = m[0]*m[4]/m[2]**2
-    rms_slopes = np.sqrt(np.pi * m[2] / 2)
-
     # Output setup
     realizations = 100
     hrms_out = np.zeros(realizations)
     alpha_out = np.zeros(realizations)
     rms_slopes_out = {"spectral": np.zeros(realizations),
                       "geometric": np.zeros(realizations)}
 
+    iso = tm.Isopowerlaw2D()
+    iso.getQ0().assign(k_l)
+    iso.getQ1().assign(k_r)
+    iso.getQ2().assign(k_s)
+    iso.getHurst().assign(Hurst)
+    sg = tm.SurfaceGeneratorFilter2D()
+    sg.setFilter(iso)
+    sg.setSizes([N, N])
+
+    # Analytical moments
+    m = iso.moments()
+    alpha = iso.alpha()
+    rms_slopes = iso.rmsSlopes()
+
     for i in range(realizations):
         # Surface generator setup
-        sg = tm.SurfaceGeneratorFilterFFT()
-        sg.getGridSize().assign(N)
-        sg.getHurst().assign(Hurst)
-        sg.getRMS().assign(hrms)
-        sg.getQ0().assign(k_l)
-        sg.getQ1().assign(k_r)
-        sg.getQ2().assign(k_s)
         sg.getRandomSeed().assign(i)
-        sg.Init()
         # Generating surface
-        surface = sg.buildSurface() / sg.analyticalRMS()
+        surface = tm.convertGrid(sg.buildSurface()) / iso.rmsHeights()
 
         # Computing RMS heights
         # We know surface.mean() = 0, so we get an unbiased estimator here
         hrms_out[i] = np.sqrt(1 / (N**2-1) * np.sum(surface**2))
 
         # Computing spectrum bandwidth
         surface /= box_size**2  # Normalize so that the computed PSD does not contain the box size
         moments = tm.SurfaceStatistics.computeMoments(surface)
         m0 = hrms_out[i]**2 / box_size**4
         m2 = moments[0]
         m4 = moments[1]
         alpha_out[i] = m0*m4/m2**2
 
         # Computing RMS slopes
         rms_slopes_out["spectral"][i] = tm.SurfaceStatistics.computeSpectralRMSSlope(surface)
         rms_slopes_out["geometric"][i] = tm.SurfaceStatistics.computeRMSSlope(surface)
 
     # print(rms_slopes_out)
     # print(rms_slopes)
     # Computing Mean Square Error
     print("""========= Analytical values =========
 hrms = {}
 alpha = {}
 m0, m2, m4 = {}
 rms_slopes = {}""".format(hrms, alpha, m, rms_slopes))
     print("========= MSE Calculations  =========")
 
     data_association = [
         (hrms_out, hrms, "hrms"),
         (alpha_out, alpha, "alpha"),
         (rms_slopes_out["spectral"], rms_slopes, "RMS Slopes Spectral"),
         (rms_slopes_out["geometric"], rms_slopes, "RMS Slopes Geometric")
     ]
 
     # Computing and printing MSE values
     mse_vals = []
     for out, analytical, string in data_association:
         mse = np.sqrt(compute_mse(out, analytical)) / analytical
         print("Normalized error (MSE) for {}: {}".format(string, mse))
         mse_vals.append(mse)
 
     # Checking tolerances
     tols = [5e-3, 5e-3, 1e-1, 1e-1]
     for mse, tol in zip(mse_vals, tols):
         if mse >= tol:
             print("{} >= {}".format(mse, tol))
             return 1
 
     tm.finalize()
     return 0
 
 
 if __name__ == "__main__":
     sys.exit(main())