diff --git a/src/core/array.hh b/src/core/array.hh
index 0fe12c0..d123773 100644
--- a/src/core/array.hh
+++ b/src/core/array.hh
@@ -1,169 +1,169 @@
 /**
  *
  * @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 "tamaas.hh"
 #include <utility>
 #include <memory>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 template <typename T>
 class Array {
 public:
   /// Default
   Array():
     _data(nullptr),
     _size(0),
     wrapped(false),
     _alloc(DEFAULT_ALLOCATOR())
   {}
   /// Empty array of given size
   Array(UInt size):
     Array(DEFAULT_ALLOCATOR()) {
     this->resize(size);
   }
   /// Copy constructor (deep)
   Array(const Array & v):
     Array() {
     this->resize(v._size);
     std::copy(v._data, v._data+v._size, this->_data);
   }
   /// Move constructor (transfers data ownership)
   Array(Array&& v):
     Array() {
-    _data = exchange(v._data, nullptr);
-    _size = exchange(v._size, 0);
-    wrapped = exchange(v.wrapped, false);
+    _data = std::exchange(v._data, nullptr);
+    _size = std::exchange(v._size, 0);
+    wrapped = std::exchange(v.wrapped, false);
   }
   /// Wrap array on data
   Array(T * data, UInt size):
     Array(DEFAULT_ALLOCATOR()) {
     this->wrapMemory(data, size);
   }
   /// Destructor
   virtual ~Array(){
     if (wrapped == false)
       _alloc.deallocate(_data, _size);
   }
 
 public:
   /// Copy operator
   Array & operator=(const Array & v){
     this->wrapped = false;
     if (v.size() != this->size()) this->resize(v.size());
 #pragma omp parallel for
     for (UInt i = 0 ; i < this->_size ; ++i) {
       _data[i] = v._data[i];
     }
     return *this;
   }
 
   /// Move operator
   Array & operator=(Array && v) {
     if (this != &v) {
       if (!wrapped) _alloc.deallocate(_data, _size);
       _data = v._data;
       _size = v._size;
       wrapped = v.wrapped;
       v._data = nullptr;
       v._size = 0;
       v.wrapped = false;
     }
     return *this;
   }
 
   /// Wrap a memory pointer
   void wrapMemory(T * data, UInt size) {
     this->_data = data;
     this->_size = size;
     this->wrapped = true;
   }
 
   /// Assign a sigle value to the array
   void assign(UInt size, const T & value){
     this->resize(size);
 #pragma omp parallel for
     for (UInt i = 0; i < size; i++) {
       _data[i] = value;
     }
   }
 
   void setWrapped(bool w) { wrapped = w; }
 
   /// Data pointer access (const)
   const T * data() const{ return _data; }
   /// Data pointer access (non-const)
   T * data() { return _data; }
   /// Data pointer setter
   void setData(T * new_ptr) { _data = new_ptr;}
   /// Resize array
   void resize(UInt size){
     if (wrapped == true)
       TAMAAS_EXCEPTION("cannot resize wrapped array");
 
     // Erase array
     if (size == 0) {
       _alloc.deallocate(_data, _size);
       this->_data = nullptr;
       this->_size = 0;
       return;
     }
 
     // Do nothing
     if (size == this->_size) return;
 
     // Allocate new data
     _alloc.deallocate(_data, _size);
 
     this->_data = _alloc.allocate(size);
     this->_size = size;
   }
   /// Access operator
   inline T & operator[] (UInt i) {
     TAMAAS_ASSERT(i < _size, "memory overflow");
     return _data[i];
   }
   /// Access operator (const)
   inline const T & operator[] (UInt i) const {
     TAMAAS_ASSERT(i < _size, "memory overflow");
     return _data[i];
   }
   /// Get size of array
   inline UInt size() const {return _size;}
 
 private:
   T * _data = nullptr;
   UInt _size = 0;
   bool wrapped = false;
   std::allocator<T> _alloc;
 };
 
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 #endif /* __ARRAY__HH__ */
diff --git a/src/core/fft_plan_manager.hh b/src/core/fft_plan_manager.hh
index 2b7060c..235a6e7 100644
--- a/src/core/fft_plan_manager.hh
+++ b/src/core/fft_plan_manager.hh
@@ -1,132 +1,132 @@
 /**
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef FFT_PLAN_MANAGER_H
 #define FFT_PLAN_MANAGER_H
 /* -------------------------------------------------------------------------- */
 #include <map>
 #include "fftransform.hh"
 #include "fftransform_fftw.hh"
 #ifdef USE_CUDA
 #include "fftransform_cufft.hh"
 #endif
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 class FFTPlanManager {
   /* ------------------------------------------------------------------------ */
   /* Constructors/Destructors                                                 */
   /* ------------------------------------------------------------------------ */
 protected:
   FFTPlanManager();
 
 public:
   virtual ~FFTPlanManager();
 
 public:
   /* ------------------------------------------------------------------------ */
   /* Methods                                                                  */
   /* ------------------------------------------------------------------------ */
 
   /// Get singleton instance
   static FFTPlanManager & get();
 
   /// Create/retrieve a plan from two surfaces
   template <UInt dim>
   FFTransform<Real, dim> & createPlan(Grid<Real, dim> & input, GridHermitian<Real, dim> & output);
 
   /// Remove all plans
   void clean();
 
   /// Destroy a plan from two surfaces
   template <UInt dim>
   void destroyPlan(Grid<Real, dim> & input, GridHermitian<Real, dim> & output);
 
   /// Destroy any plan containing a given data pointer
   template<typename T>
   void destroyPlan(const T * some);
 
   /* ------------------------------------------------------------------------ */
   /* Class Members                                                            */
   /* ------------------------------------------------------------------------ */
 private:
-  typedef std::map<std::pair<Real*,Complex*>,
-		   std::tuple<FFTransform<Real, 1> *,
-			      FFTransform<Real, 2> *,
-			      FFTransform<Real, 3> *>
-		   > FFTMap;
+  using FFTMap = std::map<std::pair<Real*,Complex*>,
+			  std::tuple<FFTransform<Real, 1> *,
+				     FFTransform<Real, 2> *,
+				     FFTransform<Real, 3> *>
+			  >;
   FFTMap plans;
   static FFTPlanManager * singleton;
 };
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 FFTransform<Real, dim> & FFTPlanManager::createPlan(Grid<Real, dim> & input,
 						    GridHermitian<Real, dim> & output){
   static_assert(dim <= 3, "Cannot do FFT of dimension higher than 3");
 
   auto index = std::make_pair(const_cast<Real*>(input.getInternalData()),
                               const_cast<Complex*>(output.getInternalData()));
 
   auto it = plans.find(index);
   auto end = plans.end();
 
   if (it == end) {
 #ifdef USE_CUDA
-    //std::get<dim-1>(plans[index]) = new FFTransformCUFFT<Real, dim>(input,output);
+    std::get<dim-1>(plans[index]) = new FFTransformFFTW<Real, dim>(input,output); // TODO cuda
 #else
     std::get<dim-1>(plans[index]) = new FFTransformFFTW<Real, dim>(input,output);
 #endif
   }
 
   return *std::get<dim-1>(plans[index]);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template<UInt dim>
 void FFTPlanManager::destroyPlan(Grid<Real, dim> & input,
                                  GridHermitian<Real, dim> & output){
 
   auto index = std::make_pair(const_cast<Real*>(input.getInternalData()),
 			      const_cast<Complex*>(output.getInternalData()));
 
   auto it = plans.find(index);
   auto end = plans.end();
 
   if (it != end){
     delete std::get<dim>(plans[index]);
     plans.erase(index);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 __END_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 #endif /* FFT_PLAN_MANAGER_H */
diff --git a/src/core/grid_base.hh b/src/core/grid_base.hh
index 373bccf..7cd0e48 100644
--- a/src/core/grid_base.hh
+++ b/src/core/grid_base.hh
@@ -1,318 +1,326 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __GRID_BASE_HH__
 #define __GRID_BASE_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "array.hh"
 #include "iterator.hh"
+#include "loop.hh"
 #include <cstddef>
 #include <utility>
 #include <vector>
 #include <limits>
 /* -------------------------------------------------------------------------- */
 
 __BEGIN_TAMAAS__
 
 /// Namespace for view index classes
 namespace view {
   /// Index class
   struct index {
     index(Int i = -1):i(i){}
     Int i;
   };
 
   /// Blocked index class
   struct blocked : public index {
     blocked(Int i):index(i){}
   };
 
   /// Free index class
   struct free : public index {
     free():index(-1){}
   };
 }
 
 template <typename T>
 class GridBase {
 public:
   /// Constructor by default
   GridBase() = default;
   /// Copy constructor
   GridBase(const GridBase & o): data(o.data),
 				nb_components(o.nb_components),
 				offset(o.offset) {}
   /// Move constructor (transfers data ownership)
   GridBase(GridBase&& o): data(std::move(o.data)),
-			  nb_components(exchange(o.nb_components, 1)),
-			  offset(exchange(o.offset, 0)) {}
+			  nb_components(std::exchange(o.nb_components, 1)),
+			  offset(std::exchange(o.offset, 0)) {}
 
   /// Destructor
   virtual ~GridBase() = default;
 
 /* -------------------------------------------------------------------------- */
 /* Iterator class                                                             */
 /* -------------------------------------------------------------------------- */
 public:
   using iterator = iterator_::iterator<T>;
   using const_iterator = iterator_::iterator<const T>;
 
 public:
   /// Compute size
   virtual UInt computeSize() const = 0;
   /// Get internal data pointer (const)
   const T * getInternalData() const {return this->data.data();}
   /// Get internal data pointer (non-const)
   T * getInternalData() {return this->data.data();}
   /// Get number of components
   UInt getNbComponents() const {return nb_components;}
   /// Set number of components
   void setNbComponents(UInt n) {nb_components = n;}
   /// Get offset
   UInt getOffset() const {return offset;}
   /// Get total size
   virtual UInt dataSize() const {return this->data.size(); }
   /// Get number of points
   UInt getNbPoints() const {return this->computeSize() / this->getNbComponents(); }
   /// Set components
   void uniformSetComponents(const GridBase<T> & vec);
 
 /* -------------------------------------------------------------------------- */
 /* Iterator methods                                                           */
 /* -------------------------------------------------------------------------- */
   virtual iterator begin(UInt n = 1) = 0;
   virtual iterator end(UInt n = 1) = 0;
   virtual const_iterator begin(UInt n = 1) const = 0;
   virtual const_iterator end(UInt n = 1) const = 0;
 
 /* -------------------------------------------------------------------------- */
 /* Operators                                                                  */
 /* -------------------------------------------------------------------------- */
+  inline T & operator()(UInt i) {
+    return this->data[i];
+  }
+
+  inline const T & operator()(UInt i) const {
+    return this->data[i];
+  }
+
 #define VEC_OPERATOR(op)					\
   template <typename T1>					\
   GridBase & operator op(const GridBase<T1> & other)		\
 
   VEC_OPERATOR(+=);
   VEC_OPERATOR(*=);
   VEC_OPERATOR(-=);
   VEC_OPERATOR(/=);
 #undef VEC_OPERATOR
 
   template <typename T1>
   T dot(const GridBase<T1> & other);
 
 #define SCALAR_OPERATOR(op)				\
   GridBase & operator op(const T & e)			\
 
   SCALAR_OPERATOR(+=);
   SCALAR_OPERATOR(*=);
   SCALAR_OPERATOR(-=);
   SCALAR_OPERATOR(/=);
   SCALAR_OPERATOR(=);
 #undef SCALAR_OPERATOR
 
   GridBase & operator=(const GridBase & o);
 
 /* -------------------------------------------------------------------------- */
 /* Min/max                                                                    */
 /* -------------------------------------------------------------------------- */
   T min() const;
   T max() const;
   T sum() const;
   T mean() const { return sum() / dataSize(); }
   T var() const;
 
 /* -------------------------------------------------------------------------- */
 /* Move/Copy                                                                  */
 /* -------------------------------------------------------------------------- */
 
 public:
   template <typename T1>
   void copy(const GridBase<T1> & other) {
     data = other.data;
     nb_components = other.nb_components;
     offset = other.offset;
   }
 
   template <typename T1>
   void move(GridBase<T1> && other) {
     data = std::move(other.data);
-    nb_components = exchange(other.nb_components, 1);
-    offset = exchange(other.offset, 0);
+    nb_components = std::exchange(other.nb_components, 1);
+    offset = std::exchange(other.offset, 0);
   }
 
 protected:
   Array<T> data;
   UInt nb_components = 1;
   UInt offset = 0;
 };
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 void GridBase<T>::uniformSetComponents(const GridBase<T> &vec) {
   if (!(vec.dataSize() == this->nb_components))
       TAMAAS_EXCEPTION("Cannot set grid field with values of vector");
 
   auto begin_it(begin()); auto end_it(end()); // silencing nvcc warnings
 
 #pragma omp parallel for
   for (auto it = begin() ; it < end_it ; ++it) {
     UInt i = it - begin_it;
     *it = vec.data[i % this->nb_components];
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 #define SCALAR_OPERATOR_IMPL(op)					\
   template <typename T>							\
   inline GridBase<T> & GridBase<T>::operator op(const T & e) {		\
-    _Pragma("omp parallel for")						\
-    for (auto it = this->begin() ; it < this->end() ; ++it) {		\
-      *it op e;								\
-    }									\
+    Loop::get().loop([&](T & val) {					\
+	val op e;							\
+      }, *this);							\
     return *this;							\
   }									\
 
 SCALAR_OPERATOR_IMPL(+=);
 SCALAR_OPERATOR_IMPL(*=);
 SCALAR_OPERATOR_IMPL(-=);
 SCALAR_OPERATOR_IMPL(/=);
 SCALAR_OPERATOR_IMPL(=);
 
 #undef SCALAR_OPERATOR_IMPL
 
 /* -------------------------------------------------------------------------- */
 
 #define VEC_OPERATOR_IMPL(op)						\
   template <typename T>							\
   template <typename T1>						\
   inline GridBase<T> & GridBase<T>::operator op(const GridBase<T1> & other) { \
     TAMAAS_ASSERT(other.computeSize() == this->computeSize(),		\
 		  "surface size does not match");			\
     auto begin_it(begin()); auto end_it(end());				\
     auto other_begin(other.begin()); auto other_it(other.begin());	\
     _Pragma("omp parallel for firstprivate(other_it)")			\
     for (auto it = begin() ; it < end_it ; ++it) {			\
       other_it = other_begin;						\
       other_it += it - begin_it;					\
       *it op *other_it;							\
     }									\
     return *this;							\
   }									\
 
 VEC_OPERATOR_IMPL(+=);
 VEC_OPERATOR_IMPL(-=);
 VEC_OPERATOR_IMPL(*=);
 VEC_OPERATOR_IMPL(/=);
 #undef VEC_OPERATOR
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 GridBase<T> & GridBase<T>::operator=(const GridBase & o) {
   this->copy(o);
   return *this;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::min() const {
   T val = std::numeric_limits<T>::max();
   auto begin_it(begin()); auto end_it(end());
 #pragma omp parallel for reduction(min:val)
   for (auto it = begin_it ; it < end_it ; ++it) {
     val = std::min(*it, val);
   }
   return val;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::max() const {
   T val = std::numeric_limits<T>::lowest();
   auto begin_it = begin(); auto end_it = end();
 #pragma omp parallel for reduction(max:val)
   for (auto it = begin_it ; it < end_it ; ++it) {
     val = std::max(*it, val);
   }
   return val;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::sum() const {
   T val = 0;
   auto begin_it(begin()); auto end_it(end());
 #pragma omp parallel for reduction(+:val)
   for (auto it = begin_it ; it < end_it ; ++it) {
     val += *it;
   }
   return val;
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 inline T GridBase<T>::var() const {
   T var = 0;
   const T mu = mean();
   auto begin_it(begin()); auto end_it(end());
 #pragma omp parallel for reduction(+:var)
   for (auto it = begin_it ; it < end_it ; ++it) {
     var += (*it - mu) * (*it - mu);
   }
 
   return var / (dataSize() - 1);
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <typename T>
 template <typename T1>
 T GridBase<T>::dot(const GridBase<T1> & other) {
   T res = 0;
   auto begin_it(begin()); auto end_it(end());
   auto other_b(other.begin()); auto other_it(other.begin());
 #pragma omp parallel for firstprivate(other_it) reduction(+: res)
   for (auto it = begin_it ; it < end_it ; ++it) {
     other_it = other_b;
     other_it += it - begin_it;
     res += *it * *other_it;
   }
 
   return res;
 }
 
 
 __END_TAMAAS__
 
 
 #endif // __GRID_BASE_HH__
diff --git a/src/core/grid_view.hh b/src/core/grid_view.hh
index 4b3ee58..7d80aaa 100644
--- a/src/core/grid_view.hh
+++ b/src/core/grid_view.hh
@@ -1,265 +1,265 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef __GRID_VIEW_HH__
 #define __GRID_VIEW_HH__
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include "grid.hh"
 #include <vector>
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 /**
  * @brief View type on grid
  * /!\ iterators do not work on this type of grid
  *
  */
 template <template <typename, UInt> class Base,
 	  typename T,
 	  UInt base_dim,
 	  UInt dim>
 class GridView : public Base<T, dim> {
 public:
   /// Constructor
   GridView(GridBase<T> & grid_base,
 	   const std::vector<view::index> & multi_index);
   /// Move constructor
   GridView(GridView&& o) :
     Base<T, dim>(std::forward<Base<T, dim>>(o)),
-    grid(exchange(o.grid, nullptr)),
+    grid(std::exchange(o.grid, nullptr)),
     indexes(std::move(o.indexes)) {}
 
   /// Destructor
   virtual ~GridView() = default;
 
 public:
   /// Compute strides and offset
   void computeStrides() override;
   /// Return view size
   UInt dataSize() const override { return this->computeSize(); }
   /// Delete resize
   void resize(const std::array<UInt, dim> &) = delete;
   /// Get a view only for the last component
   GridView lastComponentView() const;
 
 public:
   using iterator = typename GridBase<T>::iterator;
   using const_iterator = typename GridBase<T>::const_iterator;
 
   iterator begin(UInt n [[gnu::unused]] = 1) override {
     TAMAAS_EXCEPTION("iterators on views do not work as expected");
     return iterator(nullptr, 0, 0);
   }
 
   iterator end(UInt n [[gnu::unused]] = 1) override {
     TAMAAS_EXCEPTION("iterators on views do not work as expected");
     return iterator(nullptr, 0, 0);
   }
 
   const_iterator begin(UInt n [[gnu::unused]] = 1) const override {
     TAMAAS_EXCEPTION("const_iterators on views do not work as expected");
     return const_iterator(nullptr, 0, 0);
   }
 
   const_iterator end(UInt n [[gnu::unused]] = 1) const override {
     TAMAAS_EXCEPTION("const_iterators on views do not work as expected");
     return const_iterator(nullptr, 0, 0);
   }
 
 public:
   /// Access operator (1D access is extremely slow)
   template <typename... T1>
   T & operator()(T1... args);
   /// Const access operator (1D access is extremely slow)
   template <typename... T1>
   const T & operator()(T1... args) const;
 
 protected:
   Base<T, base_dim> * grid;
   std::vector<view::index> indexes;
 };
 
 /* -------------------------------------------------------------------------- */
 /* Template implementation                                                    */
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base,
 	  typename T,
 	  UInt base_dim,
 	  UInt dim>
 GridView<Base, T, base_dim, dim>::GridView(GridBase<T> & grid_base,
   const std::vector<view::index> & multi_index):
   Base<T, dim>(),
   grid(nullptr),
   indexes(multi_index)
 {
   static_assert(base_dim >= dim,
 		"view dimension must be >= than the base class");
 
   // Checking view type
   grid = dynamic_cast<Base<T, base_dim> *>(&grid_base);
 
   if (!grid) TAMAAS_EXCEPTION("given base type is incompatible with view");
 
   // Checking view integrity
   if (multi_index.size() == base_dim) {
     this->nb_components = grid->getNbComponents();
   } else if (multi_index.size() == base_dim+1) {
     this->nb_components = 1;
   } else {
     TAMAAS_EXCEPTION("Multi-index has wrong number of components");
   }
 
   UInt view_dim = 0;
   for (auto & it : multi_index) {
     if (it.i == -1) view_dim += 1;
   }
 
   if (view_dim != dim)
     TAMAAS_EXCEPTION("Multi-index has the wrong number of free components");
 
   if (indexes.size() == base_dim
       && this->nb_components == 1) indexes.push_back(view::blocked(0));
 
   this->data.wrapMemory(grid->getInternalData(), grid->dataSize());
   computeStrides();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base,
 	  typename T,
 	  UInt base_dim,
 	  UInt dim>
 void GridView<Base, T, base_dim, dim>::computeStrides() {
   /// Copying sizes with number of components in last position
   auto stride_it = this->strides.begin();
   auto size_it = this->n.begin();
   auto grid_stride_it = grid->getStrides().begin();
   auto grid_size_it = grid->sizes().begin();
 
   for (auto & index : indexes) {
     if (index.i == -1) {
       *stride_it = *grid_stride_it;
       *size_it = *grid_size_it;
       ++stride_it, ++size_it;
     } else
       this->offset += index.i * *grid_stride_it;
     ++grid_stride_it, ++grid_size_it;
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base,
 	  typename T,
 	  UInt base_dim,
 	  UInt dim>
 template <typename... T1>
 T & GridView<Base, T, base_dim, dim>::operator()(T1... args) {
   constexpr UInt nargs = sizeof...(args);
 
   static_assert(nargs == 1 || nargs == dim || nargs == dim + 1,
 		"number of arguments in operator() does not match dimension");
 
   if (nargs == 1) {
     UInt index = UInt(args...);
     std::array<UInt, dim+1> tuple = {0};
 
     tuple[dim] = index % this->nb_components;
     index -= tuple[dim];
     index /= this->nb_components;
 
     /// Computing tuple from index
     for (UInt d = dim-1 ; d > 0 ; d--) {
       tuple[d] = index % this->n[d];
       index -= tuple[d];
       index /= this->n[d];
     }
 
     tuple[0] = index;
     return Base<T, dim>::operator()(tuple);
   }
   else {
     return Base<T, dim>::operator()(args...);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base,
 	  typename T,
 	  UInt base_dim,
 	  UInt dim>
 template <typename... T1>
 const T & GridView<Base, T, base_dim, dim>::operator()(T1... args) const {
   constexpr UInt nargs = sizeof...(args);
 
   static_assert(nargs == 1 || nargs == dim || nargs == dim + 1,
 		"number of arguments in operator() does not match dimension");
 
   if (nargs == 1) {
     UInt index = UInt(args...);
     std::array<UInt, dim+1> tuple = {0};
 
     tuple[dim] = index % this->nb_components;
     index -= tuple[dim];
     index /= this->nb_components;
 
     /// Computing tuple from index
     for (UInt d = dim-1 ; d > 0 ; d--) {
       tuple[d] = index % this->n[d];
       index -= tuple[d];
       index /= this->n[d];
     }
 
     tuple[0] = index;
     return Base<T, dim>::operator()(tuple);
   }
   else {
     return Base<T, dim>::operator()(args...);
   }
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <template <typename, UInt> class Base,
 	  typename T,
 	  UInt base_dim,
 	  UInt dim>
 GridView<Base, T, base_dim, dim>
 GridView<Base, T, base_dim, dim>::lastComponentView() const {
   GridView<Base, T, base_dim, dim> tmp(*grid, indexes);
   tmp.indexes.back() = view::blocked(this->nb_components-1);
   tmp.computeStrides();
   return tmp;
 }
 
 
 __END_TAMAAS__
 
 #endif // __GRID_VIEW_HH__
diff --git a/src/core/tamaas.hh b/src/core/tamaas.hh
index 838452b..a89ee2a 100644
--- a/src/core/tamaas.hh
+++ b/src/core/tamaas.hh
@@ -1,138 +1,145 @@
 /**
  * @mainpage Welcome to Tamaas !
  *
  * @section Introduction
  * Tamaas is a spectral-boundary-element based contact library. It is made with
  * love to be fast and friendly!
  *
  * @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  * @author Valentine Rey <valentine.rey@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016-2017 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>
 #include <map>
 #include <memory>
 /* -------------------------------------------------------------------------- */
 
 #ifdef USE_CUDA
 #include "unified_allocator.hh"
 #endif
 
 /// Namespace macros
 #define __BEGIN_TAMAAS__ namespace tamaas {
 #define __END_TAMAAS__ }
 
+/// Standard definitions pulled from C++14
+namespace std {
+#if __cplusplus < 201402L
+  /// exchange
+  template<class T, class U = T>
+  T exchange(T& obj, U&& new_value)
+  {
+    T old_value = move(obj);
+    obj = forward<U>(new_value);
+    return old_value;
+  }
+  /// make_unique
+  template<typename T, typename... Args>
+  unique_ptr<T> make_unique(Args&&... args)
+  {
+    return unique_ptr<T>(new T(forward<Args>(args)...));
+  }
+#endif
+}
+
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 /// Common types definitions
 typedef double Real;
 typedef unsigned int UInt;
 typedef int Int;
 typedef std::complex<Real> Complex;
 static constexpr Real zero_threshold = 1e-14;
 
 /* -------------------------------------------------------------------------- */
 /// initialize tamaas (0 threads => let OMP_NUM_THREADS decide)
 void initialize(UInt num_threads = 0);
 /// cleanup tamaas
 void finalize();
-/// exchange implementation when C++14 is not supported
-template<class T, class U = T>
-T exchange(T& obj, U&& new_value)
-{
-    T old_value = std::move(obj);
-    obj = std::forward<U>(new_value);
-    return old_value;
-}
-/// make_unique implementation when C++14 not supported
-template<typename T, typename... Args>
-std::unique_ptr<T> make_unique(Args&&... args)
-{
-    return std::unique_ptr<T>(new T(std::forward<Args>(args)...));
-}
+
 /* -------------------------------------------------------------------------- */
 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;
 };
 /* -------------------------------------------------------------------------- */
 /// Enumeration of reduction operations
 enum operation { plus, times, min, max };
 /* -------------------------------------------------------------------------- */
 __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(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; }
 
 #ifdef USE_CUDA
 #define DEFAULT_ALLOCATOR UnifiedAllocator<T>
 #else
 #define DEFAULT_ALLOCATOR std::allocator<T>
 #endif
 
 /* -------------------------------------------------------------------------- */
 
 #endif // TAMAAS_HH
diff --git a/src/model/model_template.cpp b/src/model/model_template.cpp
index 479d666..0071c4f 100644
--- a/src/model/model_template.cpp
+++ b/src/model/model_template.cpp
@@ -1,79 +1,79 @@
 /**
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2016-2017 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 "model_template.hh"
 #include "westergaard.hh"
 /* -------------------------------------------------------------------------- */
 __BEGIN_TAMAAS__
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 ModelTemplate<type>::ModelTemplate(const std::vector<Real> & system_size,
 				   const std::vector<UInt> & discretization):
   Model(system_size, discretization),
   boundary_displacement(nullptr)
 {
   constexpr UInt dim = trait::dimension;
   constexpr UInt dim_b = trait::boundary_dimension;
   constexpr UInt nb_components = trait::components;
 
   if (system_size.size() != dim)
     TAMAAS_EXCEPTION("System size does not match model type");
   if (discretization.size() != dim)
     TAMAAS_EXCEPTION("Discretization size does not match model type");
 
   // Copying sizes for traction grid
   std::array<UInt, dim_b> traction_size;
   auto disc_it = discretization.begin();
   if (dim > dim_b) ++disc_it;
   std::copy(disc_it, discretization.end(),
 	    traction_size.begin());
 
   // Allocating
-  traction = ::tamaas::make_unique<Grid<Real, dim_b>>(traction_size, nb_components);
-  displacement = ::tamaas::make_unique<Grid<Real, dim>>(discretization, nb_components);
-  boundary_displacement = ::tamaas::make_unique<ViewType>(*displacement, trait::indices);
+  traction = std::make_unique<Grid<Real, dim_b>>(traction_size, nb_components);
+  displacement = std::make_unique<Grid<Real, dim>>(discretization, nb_components);
+  boundary_displacement = std::make_unique<ViewType>(*displacement, trait::indices);
 
-  engine = ::tamaas::make_unique<Westergaard<type>>(this);
+  engine = std::make_unique<Westergaard<type>>(this);
   engine->initInfluence();
 }
 
 /* -------------------------------------------------------------------------- */
 
 template <model_type type>
 GridBase<Real> & ModelTemplate<type>::getBoundaryDisplacement() {
   return *boundary_displacement;
 }
 
 /* -------------------------------------------------------------------------- */
 /* Template instanciation                                                     */
 /* -------------------------------------------------------------------------- */
 template class ModelTemplate<model_type::basic_1d>;
 template class ModelTemplate<model_type::basic_2d>;
 template class ModelTemplate<model_type::surface_1d>;
 template class ModelTemplate<model_type::surface_2d>;
 template class ModelTemplate<model_type::volume_1d>;
 template class ModelTemplate<model_type::volume_2d>;
 
 __END_TAMAAS__
diff --git a/tests/test_autocorrelation.py b/tests/test_autocorrelation.py
index c371133..f6a7fba 100644
--- a/tests/test_autocorrelation.py
+++ b/tests/test_autocorrelation.py
@@ -1,118 +1,120 @@
 #!/usr/bin/python
 # -*- coding: utf-8 -*-
 #
 # @author Valentine Rey <valentine.rey@epfl.ch>
 #
 # @section LICENSE
 #
 # Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
 # Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
 # Solides)
 #
 # Tamaas is free  software: you can redistribute it and/or  modify it under the
 # terms  of the  GNU Lesser  General Public  License as  published by  the Free
 # Software Foundation, either version 3 of the License, or (at your option) any
 # later version.
 #
 # Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 # A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
 # details.
 #
 # You should  have received  a copy  of the GNU  Lesser General  Public License
 # along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
 #
 #
 ################################################################
 from __future__ import print_function
 
 from tamaas import *
 
 import numpy as np
 from scipy.special import jv
 import scipy.integrate as integrate
 from numpy.fft import *
 
 import sys, imp
 
 
 def spectrum(r, kl, kr, ks, C, H):
     if r > ks:
         return 0
     elif r > kr:
         return C * (r/float(kr))**(-2.*(H+1))
     elif r > kl:
         return C
     else:
         return 0
 
 
 # n is normalization parameter
 def integrate_bessel(x, spectrum, n):
     result = 2 * np.pi * integrate.quad(lambda y: y*jv(0, y*x) * spectrum(n*y), 0, np.inf)[0]
     return result
 
 
 def main():
+    initialize()
     grid_size = 256
 
     lambda_l = 8.
     lambda_r = 8.
     lambda_s = 1.
 
     L = 4 * lambda_l
 
     kl = int(L / lambda_l)
     kr = int(L / lambda_r)
     ks = int(L / lambda_s)
 
 
 
     # Computing ACF from integral with Bessel function
     l = grid_size/2 # points on curve
     x = np.linspace(0, L/2, l)
     bessel_acf = np.zeros(l)
     integrate_bessel_v = np.vectorize(integrate_bessel, otypes=[np.float])
     bessel_acf = (L / (2*np.pi) )**2 *integrate_bessel_v(x, lambda x: spectrum(x,
                                                                 kl, kr, ks,
                                                                 1., 0.8), L / (2*np.pi))
 
 
     nsamples = 1000
     etages=np.zeros((grid_size,grid_size,nsamples))
 
     # Computing ACF with Hu & Tonder algo reimplemented
 
     for i in range(nsamples):
             #generate surface
         SG = SurfaceGeneratorFilterFFT()
         SG.getGridSize().assign(grid_size)
         SG.getHurst().assign(0.8)
         SG.getQ0().assign(kl);
         SG.getQ1().assign(kr);
         SG.getQ2().assign(ks);
         SG.getRandomSeed().assign(i);
         SG.Init()
         a = SG.buildSurface()
         etages[:,:,i]=a
 
 
 
     cov=np.zeros(grid_size/2)
 
     for i in np.arange(grid_size/2,grid_size,1):
         cov[i-grid_size/2]=np.sum((etages[grid_size/2,i,:]-np.mean(etages[grid_size/2,i,:]))*(etages[grid_size/2,grid_size/2,:]-np.mean(etages[grid_size/2,grid_size/2,:])))/nsamples
 
 
 
     error=np.sum((cov-bessel_acf)*(cov-bessel_acf))/np.sum(bessel_acf*bessel_acf)
     print(error)
 
     if error > 1e-1 :
             return 1
 
 
+    finalize()
     return 0
 
 if __name__ == "__main__":
     sys.exit(main())
diff --git a/tests/test_fftransform.py b/tests/test_fftransform.py
index e678cce..7c0289b 100644
--- a/tests/test_fftransform.py
+++ b/tests/test_fftransform.py
@@ -1,100 +1,102 @@
 #!/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
 
 import sys
 import tamaas as tm
 import numpy as np
 from numpy.linalg import norm
 
 
 def main():
+    tm.initialize()
     size = 128
     surface = np.random.randn(128, 128)
 
     # Compute fourier transform through numpy
     numpy_spectral = np.fft.rfft2(surface)
 
     # Compute fourier tranform through tamaas
     tamaas_spectral = np.zeros((size, size/2+1), np.complex)
     transform = tm.FFTransformFFTWReal(surface, tamaas_spectral)
     transform.forwardTransform()
 
     spectral = tamaas_spectral
     error_npy = np.abs(np.sum((spectral-numpy_spectral)**2))
 
     surface_copy = surface.copy()
     surface[:, :] = 0
 
     transform.backwardTransform()
     error_back = np.abs(np.sum((surface_copy-surface)**2))
 
     print(error_back)
     print(error_npy)
 
     if error_back > 1e-16 or error_npy > 1e-16:
         return 1
 
     # Testing multi-component FFT
     multi = np.zeros((size, size, 3))
     hermit = np.zeros((size, size/2+1, 3), np.complex)
     transform = tm.FFTransformFFTWReal(multi, hermit)
     norms = []
 
     # Checking no overflow in other components of spectral representation
     multi.fill(0.)
     multi[:, :, 0] = surface[:, :]
     transform.forwardTransform()
     norms.append(norm(hermit[:, :, (1, 2)]))
 
     multi.fill(0.)
     multi[:, :, 1] = surface[:, :]
     transform.forwardTransform()
     norms.append(norm(hermit[:, :, (0, 2)]))
 
     multi.fill(0.)
     multi[:, :, 2] = surface[:, :]
     transform.forwardTransform()
     norms.append(norm(hermit[:, :, (0, 1)]))
 
     for i in range(3):
         if norms[i] != 0:
             print("Error in multi-component data separation")
             return 1
         multi[:, :, i] = surface[:, :]
 
     # Multi is full of surface
     transform.forwardTransform()
     for i in range(3):
         error = norm(hermit[:, :, i] - numpy_spectral) / norm(numpy_spectral)
         if error > 1e-10:
             print("Error in multi-component spectral data : {} => {}".format(i, error))
             return 1
 
+    tm.finalize()
     return 0
 
 
 if __name__ == '__main__':
     sys.exit(main())
diff --git a/tests/test_flood_fill.py b/tests/test_flood_fill.py
index e1929c3..aeaeaf1 100644
--- a/tests/test_flood_fill.py
+++ b/tests/test_flood_fill.py
@@ -1,69 +1,73 @@
 #!/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 tamaas as tm
 import numpy as np
 
+tm.initialize()
+
 field = np.zeros((18, 18))
 
 # Cluster of permeter 18 area 13
 field[2:4, 3:6] = 1.
 field[4  , 4:6] = 1.
 field[5  , 5  ] = 1.
 field[3:5, 6:8] = 1.
 
 # Cluster of perimeter 8 area 4
 field[14:16, 3:5] = 1.
 
 # Cluster of perimeter 20 area 11
 field[9:11, 8:11 ] = 1.
 field[7:9 , 9    ] = 1.
 field[10  , 11:14] = 1.
 
 # Cluster of perimeter 18 area 9
 field[3:5, 14:16] = 1.
 field[5:10, 15] = 1.
 
 bool_field = np.zeros_like(field, np.bool)
 bool_field[field > 0] = True
 
 clusters = tm.FloodFill.getClusters(bool_field, False)
 
 
 # Dummy class for clusters
 class dummy:
     def __init__(self, area, perimeter):
         self.area = area
         self.perimeter = perimeter
 
 
 # Solution
 ref = [dummy(13, 18), dummy(9, 18), dummy(11, 20), dummy(4, 8)]
 
 for x, y in zip(clusters, ref):
     if x.getArea() != y.area:
         sys.exit(1)
+
+tm.finalize()
diff --git a/tests/test_hertz_adhesion.py b/tests/test_hertz_adhesion.py
index 9f6fa49..10d181f 100644
--- a/tests/test_hertz_adhesion.py
+++ b/tests/test_hertz_adhesion.py
@@ -1,129 +1,131 @@
 # -*- coding: utf-8 -*-
 # -----------------------------------------------------------------------------
 # @author Valentine Rey <valentine.rey@epfl.ch>
 #
 # @section LICENSE
 #
 # Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
 # Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
 # Solides)
 #
 # Tamaas is free  software: you can redistribute it and/or  modify it under the
 # terms  of the  GNU Lesser  General Public  License as  published by  the Free
 # Software Foundation, either version 3 of the License, or (at your option) any
 # later version.
 #
 # Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 # A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
 # details.
 #
 # You should  have received  a copy  of the GNU  Lesser General  Public License
 # along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
 # -----------------------------------------------------------------------------
 from __future__ import print_function, division
 
 import sys
 import tamaas as tm
 import numpy as np
 import scipy.optimize as scio
 
 
 def constructHertzProfile(size, curvature):
     radius = 1. / curvature
     x = np.linspace(-0.5, 0.5, size)
     y = np.linspace(-0.5, 0.5, size)
     x, y = np.meshgrid(x, y)
     surface = -(x**2 + y**2)/(2*radius)
     return surface.copy()
 
 
 def main():
+    tm.initialize()
     grid_size = 512
     curvature = 0.5
     effective_modulus = 1.
     surface_energy = 0.0001
     rho = 0.002
     radius = 1./curvature
     K = 4*effective_modulus/3
     KI = np.sqrt(surface_energy*2*effective_modulus)
     tabor = surface_energy**(2/3)/(rho*effective_modulus**(2/3))*(1/curvature)**(1/3)
     sig0 = surface_energy/(rho)
     lam = 2 * sig0/(np.pi*surface_energy*K**2/radius)**(1/3)
 
     surface = constructHertzProfile(grid_size, curvature)
     bem = tm.BemGigipol(surface)
     bem.setDumpFreq(100)
     functional = tm.MaugisAdhesionFunctional(bem)
     functional.setParameter('rho', rho)
     functional.setParameter('surface_energy', surface_energy)
     bem.setEffectiveModulus(effective_modulus)
     bem.addFunctional(functional)
     bem.computeEquilibrium(1e-8, 0.00116480372339)
 
     tractions = bem.getTractions()
     true_gap = bem.getGap()
     true_displacements = bem.getTrueDisplacements()
 
     adhesive_pressure = np.zeros((grid_size, grid_size))
     for i in np.arange(0, grid_size, 1):
         for j in np.arange(0, grid_size, 1):
             if (rho > true_gap[i, j]):
                 adhesive_pressure[i, j] = surface_energy/rho
 
     # computed solution
     p0 = np.mean(tractions)
     p_max = np.max(tractions)
     contact_points = np.where(true_gap == 0)
     contact_area = np.size(contact_points)/2/(float(grid_size)**2)
     ag = (contact_area/np.pi)**0.5
     adh_points = np.where(adhesive_pressure != 0.)
     adh_area = np.size(adh_points)/2/(float(grid_size)**2)
     cg = (adh_area/np.pi)**0.5
     delta = ag**2/radius-8*sig0/(3*K)*np.sqrt(cg**2-ag**2)
 
     # Exact solution
     beta = (1-np.exp(lam/-0.924))/1.02
     a_chap = 1.54+((2.28*lam**1.3-1)/(2.28*lam**1.3+1))*0.279
     L_chap = -7./4.+((4.04*lam**1.4-1)/(4.04*lam**1.4+1))/4.
     a0 = a_chap*(np.pi*surface_energy*radius**2/K)**(1./3.)
     pc = L_chap*np.pi*surface_energy*radius
     P = -0.00078
     ac = a0*((beta+np.sqrt(1-P/pc))/(1+beta))**(2./3.)
 
     def fonction1(m):
         return np.sqrt(m**2-1)+m**2*np.arctan(np.sqrt(m**2-1))-np.pi*KI/(np.sqrt(np.pi*ac)*sig0)
     m = scio.brenth(fonction1, 1., 2.)
     cc = m*ac
     delta_cos = ac**2/radius-8*sig0/(3*K)*np.sqrt(cc**2-ac**2)
 
     x = np.linspace(-0.5, 0.5, grid_size)
     y = np.linspace(-0.5, 0.5, grid_size)
     x, y = np.meshgrid(x, y)
     a = ac
     c = cc
     # stress
     sig_hertz = np.zeros((grid_size, grid_size))
     disp = np.zeros((grid_size, grid_size))
     for i in np.arange(0, grid_size, 1):
         for j in np.arange(0, grid_size, 1):
             if a**2>(x[i,j]**2 + y[i,j]**2):
                 disp[i,j]=3*K/(2*np.pi*radius)*np.sqrt(a**2-(x[i,j]**2 + y[i,j]**2))-2*sig0/np.pi*np.arctan(np.sqrt((c**2-a**2)/(a**2-(x[i,j]**2 + y[i,j]**2))))
             if c**2>(x[i,j]**2 + y[i,j]**2)>a**2:
                 disp[i,j]=-sig0
     sig_hertz = disp
 
     a_error = abs((ac-ag)/ac)
     c_error = abs((cc-cg)/cc)
     delta_error = abs((delta-delta_cos)/delta)
     pressure_error = abs((p0-np.mean(sig_hertz))/np.mean(sig_hertz))
 
     if a_error > 1.e-2 or c_error > 1.2e-2 or delta_error > 5.e-2 or pressure_error > 1.e-1:
         for err in [a_error, c_error, delta_error, pressure_error]:
             print(err)
         return 1
+    tm.finalize()
     return 0
 
 
 if __name__ == "__main__":
     sys.exit(main())
diff --git a/tests/test_hertz_disp.py b/tests/test_hertz_disp.py
index fd42265..6396ab9 100644
--- a/tests/test_hertz_disp.py
+++ b/tests/test_hertz_disp.py
@@ -1,136 +1,138 @@
 # -*- coding: utf-8 -*-
 # -----------------------------------------------------------------------------
 # @author Lucas Frérot <lucas.frerot@epfl.ch>
 # @author Valentine Rey <valentine.rey@epfl.ch>
 #
 # @section LICENSE
 #
 # Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
 # Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
 # Solides)
 #
 # Tamaas is free  software: you can redistribute it and/or  modify it under the
 # terms  of the  GNU Lesser  General Public  License as  published by  the Free
 # Software Foundation, either version 3 of the License, or (at your option) any
 # later version.
 #
 # Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 # A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
 # details.
 #
 # You should  have received  a copy  of the GNU  Lesser General  Public License
 # along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
 # -----------------------------------------------------------------------------
 from __future__ import print_function, division
 
 import sys
 import tamaas as tm
 
 import numpy as np
 
 
 def constructHertzProfile(size, curvature):
     radius = 1. / curvature
     x = np.linspace(-0.5, 0.5, size)
     y = np.linspace(-0.5, 0.5, size)
     x, y = np.meshgrid(x, y)
     surface = np.sqrt(radius**2 - x**2 - y**2)
     surface -= surface.min()
     return surface.copy()
 
 
 def computeHertzDisplacement(e_star, contact_size, max_pressure, size):
     x = np.linspace(-0.5, 0.5, size)
     y = np.linspace(-0.5, 0.5, size)
     x, y = np.meshgrid(x, y)
     disp = np.pi * max_pressure / (4 * contact_size * e_star) * (2 * contact_size**2 -
                                                                  (x**2 + y**2))
     return disp.copy()
 
 
 def main():
+    tm.initialize()
     grid_size = 512
     curvature = 0.5
     effective_modulus = 1.
     load = 0.12
 
     surface = constructHertzProfile(grid_size, curvature)
     model = tm.ModelFactory.createModel(tm.model_type_basic_2d,
                                         [1., 1.],
                                         [grid_size, grid_size])
     model.setElasticity(1, 0)
     solver = tm.PolonskyKeerRey(model, surface, 1e-12,
                                 tm.PolonskyKeerRey.gap,
                                 tm.PolonskyKeerRey.gap)
 
     solver.solve(load - surface.mean())
 
     tractions = model.getTraction()
     true_displacements = model.getDisplacement()
     true_gap = true_displacements - surface
     pressure = np.mean(tractions)
     contact_points = np.where(1.0e-10 >= true_gap)
     contact_area = (np.size(contact_points)/2)/float(grid_size*grid_size)
 
     print('The contact area computed with the gap map is '+str(contact_area))
 
     hertz_contact_size = (3 * pressure / (4 * curvature * effective_modulus))**(1. / 3.)
 
     print('Exact Hertz contact radius is '+str(hertz_contact_size))
     hertz_area = np.pi * hertz_contact_size**2
     print('Exact Hertz contact area is '+str(hertz_area))
     area_error = np.abs(hertz_area - contact_area) / hertz_area
     print("Area error: {}".format(area_error))
 
     # Testing maximum pressure
     max_pressure = tractions.max()
     hertz_max_pressure = (6 * pressure * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi
     pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure
     print('Exact Hertz max pressure is '+str(hertz_max_pressure))
     print('max pressure is '+str(np.max(tractions)))
     print("Max pressure error: {}".format(pressure_error))
 
     # Testing displacements
     hertz_displacements = computeHertzDisplacement(effective_modulus,
                                                    hertz_contact_size,
                                                    hertz_max_pressure,
                                                    grid_size)
 
     # Selecing only the points that are in contact
     contact_indexes = [(i, j)
                        for i in range(grid_size) for j in range(grid_size)
                        if true_gap[i, j] < 1e-12]
 
     # Displacements of bem are centered around the mean of the whole surface
     # and Hertz displacements are not centered, so we need to compute mean
     # on the contact zone for both arrays
     bem_mean = 0.
     hertz_mean = 0.
     for index in contact_indexes:
         bem_mean += true_displacements[index]
         hertz_mean += hertz_displacements[index]
 
     bem_mean /= len(contact_indexes)
     hertz_mean /= len(contact_indexes)
     # Correction applied when computing error
     correction = hertz_mean - bem_mean
 
     # Computation of error of displacement in contact zone
     error = 0.
     hertz_norm = 0.
     for index in contact_indexes:
         error += (hertz_displacements[index] - true_displacements[index] - correction)**2
         hertz_norm += (hertz_displacements[index] - hertz_mean)**2
 
     displacement_error = np.sqrt(error / hertz_norm)
     print("Displacement error (in contact zone): {}".format(displacement_error))
 
     if area_error > 1e-2 or pressure_error > 1e-2 or displacement_error > 2e-3:
         return 1
 
+    tm.finalize()
     return 0
 
 
 if __name__ == "__main__":
     sys.exit(main())
diff --git a/tests/test_hertz_kato.py b/tests/test_hertz_kato.py
index 3e14405..fb2d529 100644
--- a/tests/test_hertz_kato.py
+++ b/tests/test_hertz_kato.py
@@ -1,122 +1,124 @@
 #!/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 division, print_function
 
 import sys
 import numpy as np
 import tamaas as tm
 
 
 def constructHertzProfile(size, curvature):
     radius = 1. / curvature
     x = np.linspace(-0.5, 0.5, size)
     y = np.linspace(-0.5, 0.5, size)
     x, y = np.meshgrid(x, y)
     surface = np.sqrt(radius**2 - x**2 - y**2)
     surface -= surface.mean()
     return surface.copy()
 
 
 def computeHertzDisplacement(e_star, contact_size, max_pressure, size):
     x = np.linspace(-0.5, 0.5, size)
     y = np.linspace(-0.5, 0.5, size)
     x, y = np.meshgrid(x, y)
     disp = np.pi * max_pressure / (4 * contact_size * e_star) * (2 * contact_size**2 -
                                                                  (x**2 + y**2))
     return disp.copy()
 
 
 def main():
+    tm.initialize()
     grid_size = 512
     curvature = 0.1
     effective_modulus = 1.
     load = 0.0001
 
     surface = constructHertzProfile(grid_size, curvature)
 
     bem = tm.BemKato(surface)
     bem.setEffectiveModulus(effective_modulus)
     bem.setMaxIterations(100000)
     bem.computeEquilibrium(1e-7, load)
 
     tractions = bem.getTractions()
     displacements = bem.getDisplacements()
 
     # Testing contact area against Hertz solution for solids of revolution
     contact_area = tm.SurfaceStatistics.computeContactArea(tractions)
     hertz_contact_size = (3 * load / (4 * curvature * effective_modulus))**(1. / 3.)
     hertz_area = np.pi * hertz_contact_size**2
     area_error = np.abs(hertz_area - contact_area) / hertz_area
     print("Area error: {}".format(area_error))
 
     # Testing maximum pressure
     max_pressure = tractions.max()
     hertz_max_pressure = (6 * load * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi
     pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure
     print("Max pressure error: {}".format(pressure_error))
 
     # Testing displacements
     hertz_displacements = computeHertzDisplacement(effective_modulus,
                                                    hertz_contact_size,
                                                    hertz_max_pressure,
                                                    grid_size)
 
     # Selecing only the points that are in contact
     contact_indexes = [(i, j)
                        for i in range(grid_size) for j in range(grid_size)
                        if tractions[i, j] > 0]
 
     # Displacements of bem are centered around the mean of the whole surface
     # and Hertz displacements are not centered, so we need to compute mean
     # on the contact zone for both arrays
     bem_mean = 0.
     hertz_mean = 0.
     for index in contact_indexes:
         bem_mean += displacements[index]
         hertz_mean += hertz_displacements[index]
 
     bem_mean /= len(contact_indexes)
     hertz_mean /= len(contact_indexes)
     # Correction applied when computing error
     correction = hertz_mean - bem_mean
 
     # Computation of error of displacement in contact zone
     error = 0.
     hertz_norm = 0.
     for index in contact_indexes:
         error += (hertz_displacements[index] - displacements[index] - correction)**2
         hertz_norm += (hertz_displacements[index] - hertz_mean)**2
 
     displacement_error = np.sqrt(error / hertz_norm)
     print("Displacement error (in contact zone): {}".format(displacement_error))
 
     if area_error > 6e-3 or pressure_error > 4e-3 or displacement_error > 3e-4:
         return 1
 
+    tm.finalize()
     return 0
 
 
 if __name__ == "__main__":
     sys.exit(main())
diff --git a/tests/test_hertz_pressure.py b/tests/test_hertz_pressure.py
index cfe30bf..9d5118e 100644
--- a/tests/test_hertz_pressure.py
+++ b/tests/test_hertz_pressure.py
@@ -1,124 +1,126 @@
 #!/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 division, print_function
 
 import sys
 import numpy as np
 import tamaas as tm
 
 
 def constructHertzProfile(size, curvature):
     radius = 1. / curvature
     x = np.linspace(-0.5, 0.5, size)
     y = np.linspace(-0.5, 0.5, size)
     x, y = np.meshgrid(x, y)
     surface = np.sqrt(radius**2 - x**2 - y**2)
     surface -= surface.mean()
     return surface.copy()
 
 
 def computeHertzDisplacement(e_star, contact_size, max_pressure, size):
     x = np.linspace(-0.5, 0.5, size)
     y = np.linspace(-0.5, 0.5, size)
     x, y = np.meshgrid(x, y)
     disp = np.pi * max_pressure / (4 * contact_size * e_star) * (2 * contact_size**2 -
                                                                  (x**2 + y**2))
     return disp.copy()
 
 
 def main():
+    tm.initialize()
     grid_size = 1024
     curvature = 0.1
     effective_modulus = 1.
     load = 0.0001
 
     surface = constructHertzProfile(grid_size, curvature)
     model = tm.ModelFactory.createModel(tm.model_type_basic_2d,
                                         [1., 1.],
                                         [grid_size, grid_size])
     model.setElasticity(1, 0)
     solver = tm.PolonskyKeerRey(model, surface, 1e-12,
                                 tm.PolonskyKeerRey.pressure,
                                 tm.PolonskyKeerRey.pressure)
     solver.solve(load)
 
     tractions = model.getTraction()
     displacements = model.getDisplacement()
 
     # Testing contact area against Hertz solution for solids of revolution
     contact_area = tm.SurfaceStatistics.computeContactArea(tractions)
     hertz_contact_size = (3 * load / (4 * curvature * effective_modulus))**(1. / 3.)
     hertz_area = np.pi * hertz_contact_size**2
     area_error = np.abs(hertz_area - contact_area) / hertz_area
     print("Area error: {}".format(area_error))
 
     # Testing maximum pressure
     max_pressure = tractions.max()
     hertz_max_pressure = (6 * load * effective_modulus**2 * curvature ** 2)**(1. / 3.) / np.pi
     pressure_error = np.abs(hertz_max_pressure - max_pressure) / hertz_max_pressure
     print("Max pressure error: {}".format(pressure_error))
 
     # Testing displacements
     hertz_displacements = computeHertzDisplacement(effective_modulus,
                                                    hertz_contact_size,
                                                    hertz_max_pressure,
                                                    grid_size)
 
     # Selecing only the points that are in contact
     contact_indexes = [(i, j) for i in range(grid_size) for j in range(grid_size)
                        if tractions[i, j] > 0]
 
     # Displacements of bem are centered around the mean of the whole surface
     # and Hertz displacements are not centered, so we need to compute mean
     # on the contact zone for both arrays
     bem_mean = 0.
     hertz_mean = 0.
     for index in contact_indexes:
         bem_mean += displacements[index]
         hertz_mean += hertz_displacements[index]
 
     bem_mean /= len(contact_indexes)
     hertz_mean /= len(contact_indexes)
     # Correction applied when computing error
     correction = hertz_mean - bem_mean
 
     # Computation of error of displacement in contact zone
     error = 0.
     hertz_norm = 0.
     for index in contact_indexes:
         error += (hertz_displacements[index] - displacements[index] - correction)**2
         hertz_norm += (hertz_displacements[index] - hertz_mean)**2
 
     displacement_error = np.sqrt(error / hertz_norm)
     print("Displacement error (in contact zone): {}".format(displacement_error))
 
     if area_error > 1e-3 or pressure_error > 3e-3 or displacement_error > 1e-4:
         return 1
 
+    tm.finalize()
     return 0
 
 
 if __name__ == "__main__":
     sys.exit(main())
diff --git a/tests/test_patch_westergaard.py b/tests/test_patch_westergaard.py
index 6f69f73..4983210 100644
--- a/tests/test_patch_westergaard.py
+++ b/tests/test_patch_westergaard.py
@@ -1,74 +1,76 @@
 #!/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 division, print_function
 
 import sys
 import numpy as np
 from numpy.linalg import norm
 import tamaas as tm
 
 
 def constructSinProfile(size, mode, amplitude):
     x = np.linspace(0, 1, size, endpoint=False)
     y = np.linspace(0, 1, size, endpoint=False)
     x, y = np.meshgrid(x, y)
     surface = amplitude * np.sin(2*np.pi*x*mode) * np.sin(2*np.pi*y*mode)
     return surface.copy()
 
 
 def main():
+    tm.initialize()
     E = 3.
     nu = 0.
     E_star = E / (1-nu**2)
     grid_size = 4
     input_pressure = constructSinProfile(grid_size, 1, 1)
     model = tm.ModelFactory.createModel(tm.model_type_basic_2d,
                                         [1., 1.],
                                         [grid_size, grid_size])
     model.setElasticity(E, nu)
     model.getTraction()[:, :] = input_pressure[:, :]
     model.solveNeumann()
     output_displ = model.getDisplacement()
 
     solution = constructSinProfile(grid_size, 1, 1 / (np.pi * E_star * np.sqrt(2)))
     error = norm(solution - output_displ) / norm(solution)
 
     if error > 1e-17:
         print("Neumann error = {}".format(error))
         return 1
 
     output_displ[:, :] = solution[:, :]
     model.solveDirichlet()
     error = norm(input_pressure - model.getTraction()) / norm(input_pressure)
 
     if error > 2e-17:
         print("Dirichlet error = {}".format(error))
         return 1
+    tm.finalize()
     return 0
 
 
 if __name__ == "__main__":
     sys.exit(main())
diff --git a/tests/test_saturated_pressure.py b/tests/test_saturated_pressure.py
index ec3f91e..32fcb17 100644
--- a/tests/test_saturated_pressure.py
+++ b/tests/test_saturated_pressure.py
@@ -1,80 +1,82 @@
 #!/usr/bin/python
 # -*- coding: utf-8 -*-
 #
 # @author Valentine Rey <valentine.rey@epfl.ch>
 #
 # @section LICENSE
 #
 # Copyright (©)  2016 EPFL  (Ecole Polytechnique  Fédérale de
 # Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
 # Solides)
 #
 # Tamaas is free  software: you can redistribute it and/or  modify it under the
 # terms  of the  GNU Lesser  General Public  License as  published by  the Free
 # Software Foundation, either version 3 of the License, or (at your option) any
 # later version.
 #
 # Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
 # WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
 # A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
 # details.
 #
 # You should  have received  a copy  of the GNU  Lesser General  Public License
 # along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
 #
 #
 ################################################################
 from __future__ import print_function, division
 
 import sys, imp
 import numpy as np
 import tamaas as tm
 
 
 def main():
+    tm.initialize()
     grid_size = 256
     load = 0.06
     effective_modulus = 1.
     p_sat = 0.4
     do_plot = 1
 
     SG = tm.SurfaceGeneratorFilterFFT()
     SG.getGridSize().assign(grid_size)
     SG.getHurst().assign(0.8)
     SG.getRMS().assign(0.01)
     SG.getQ0().assign(4)
     SG.getQ1().assign(4)
     SG.getQ2().assign(16)
     SG.getRandomSeed().assign(2)
     SG.Init()
     surface = SG.buildSurface()
     surface -= np.max(surface)
 
     bem = tm.BemKato(surface)
     bem.setDumpFreq(100)
     bem.setEffectiveModulus(effective_modulus)
     bem.setMaxPressure(p_sat)
     bem.setMaxIterations(6000)
 
     bem.computeEquilibrium(1e-12,load)
     tractions = bem.getTractions()
     true_displacements = bem.getTrueDisplacements()
     true_gap = true_displacements-surface
 
     mean_pressure_error = abs(np.mean(tractions)-load)/load
     max_pressure_error = np.max(tractions)-p_sat
 
     gaps = true_gap
     gaps[np.where(tractions == p_sat)] = 0.
     press = tractions
     orthogonality_error = np.max(press*gaps) / (press.max() * gaps.max())
 
     if mean_pressure_error > 1.e-12 or max_pressure_error > 1.e-12 or orthogonality_error > 1.e-10:
         print(mean_pressure_error, max_pressure_error, orthogonality_error)
         return 1
 
+    tm.finalize()
     return 0
 
 
 if __name__ == "__main__":
     sys.exit(main())
diff --git a/tests/test_westergaard.py b/tests/test_westergaard.py
index 63b66b6..7af9a55 100644
--- a/tests/test_westergaard.py
+++ b/tests/test_westergaard.py
@@ -1,115 +1,117 @@
 #!/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 division, print_function
 
 import sys
 import numpy as np
 from numpy.linalg import norm
 from numpy import sin, cos, pi, sqrt, log
 import tamaas as tm
 
 
 def constructSinProfile(size, mode, amplitude):
     x = np.linspace(0, 1, size, endpoint=False)
     y = np.linspace(0, 1, size, endpoint=False)
     x, y = np.meshgrid(x, y)
     surface = amplitude * np.sin(2*np.pi*x*mode)
     return surface.copy()
 
 
 def main():
+    tm.initialize()
     grid_size = 2187  # Power of 3 for better accuracy
     lamda = 1.
     delta = 0.1
     effective_modulus = 1.
 
     p_star = np.pi * effective_modulus * delta / lamda  # Full contact load
 
     load = 0.1 * p_star
 
     surface = constructSinProfile(grid_size, 1 / lamda, delta)
 
     model = tm.ModelFactory.createModel(tm.model_type_basic_2d,
                                         [1., 1.],
                                         [grid_size, grid_size])
     model.setElasticity(1, 0)
     solver = tm.PolonskyKeerRey(model, surface, 1e-12,
                                 tm.PolonskyKeerRey.pressure,
                                 tm.PolonskyKeerRey.pressure)
     solver.setDumpFrequency(1)
     solver.solve(load)
 
     tractions = model.getTraction()
     displacements = model.getDisplacement()
 
     # Testing contact area against Westergaard solution
     contact_area = np.sum(tractions > 0) / grid_size**2
     westergaard_contact_size = 2 * lamda / np.pi * np.arcsin(np.sqrt(load / p_star))
     area_error = np.abs(contact_area - westergaard_contact_size / lamda) / \
         (westergaard_contact_size / lamda)
     print("Area error: {}".format(area_error))
 
     # Testing agains pressure from Westergaard solution
     x = np.linspace(0, 1, grid_size, endpoint=False)
     y = np.linspace(0, 1, grid_size, endpoint=False)
     x, y = np.meshgrid(x, y)
     x -= 1 / 4
     a = westergaard_contact_size / 2
     pressure = np.zeros_like(tractions)
     pressure = 2 * load * np.cos(np.pi*x/lamda)/np.sin(np.pi*a/lamda)**2 * \
         np.sqrt(np.sin(np.pi*a/lamda)**2 - np.sin(np.pi*x/lamda)**2)
     pressure[np.abs(x) > a] = 0
 
     pressure_error = norm(pressure - tractions) / norm(pressure)
     print("Pressure error: {}".format(pressure_error))
 
     # Energetic error (displacement reference in Johnson)
     displacements -= displacements.max()
 
     psi = np.pi * np.abs(x) / lamda
     psi_a = np.pi * a / lamda
 
     disp_w = load * lamda / (pi * effective_modulus * sin(psi_a)) * (
         cos(2*psi) + 2*sin(psi)*sqrt(sin(psi)**2 - sin(psi_a)**2)
         - 2 * sin(psi_a)**2 * log((sin(psi) + sqrt(sin(psi)**2 - sin(psi_a)**2)) / sin(psi_a)))
     disp_w[np.abs(x) < a] = load * lamda / (pi * effective_modulus * sin(psi_a)) * \
         cos(2*psi[np.abs(x) < a])
 
     disp_w -= disp_w.max()
 
     energy_error = np.sum((tractions-pressure)*(displacements-disp_w)) / \
         norm(pressure * disp_w)/grid_size**2
 
     print("Energy error: {}".format(energy_error))
 
     if area_error > 1e-4 or pressure_error > 5e-4 or energy_error > 2e-8:
         return 1
 
+    tm.finalize()
     return 0
 
 
 if __name__ == "__main__":
     sys.exit(main())