diff --git a/src/core/allocator.hh b/src/core/allocator.hh
new file mode 100644
index 0000000..fb9f783
--- /dev/null
+++ b/src/core/allocator.hh
@@ -0,0 +1,44 @@
+/*
+ *  SPDX-License-Indentifier: AGPL-3.0-or-later
+ *
+ *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
+ *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
+ *  Copyright (©) 2020-2022 Lucas Frérot
+ *
+ *  This program is free software: you can redistribute it and/or modify
+ *  it under the terms of the GNU Affero General Public License as published
+ *  by the Free Software Foundation, either version 3 of the License, or
+ *  (at your option) any later version.
+ *
+ *  This program 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 Affero General Public License for more details.
+ *
+ *  You should have received a copy of the GNU Affero General Public License
+ *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
+ *
+ */
+/* -------------------------------------------------------------------------- */
+#ifndef TAMAAS_ALLOCATOR_HH
+#define TAMAAS_ALLOCATOR_HH
+
+#include "tamaas.hh"
+
+#ifdef TAMAAS_USE_CUDA
+#include "cuda/unified_allocator.hh"
+#else
+#include "fftw/fftw_allocator.hh"
+#endif
+
+namespace tamaas {
+#ifdef TAMAAS_USE_CUDA
+template <typename T>
+using Allocator = UnifiedAllocator<T>;
+#else
+template <typename T>
+using Allocator = FFTWAllocator<T>;
+#endif
+}  // namespace tamaas
+
+#endif  // TAMAAS_ALLOCATOR_HH
diff --git a/src/core/array.hh b/src/core/array.hh
index c145d01..c4880ab 100644
--- a/src/core/array.hh
+++ b/src/core/array.hh
@@ -1,179 +1,180 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef ARRAY_HH
 #define ARRAY_HH
 /* -------------------------------------------------------------------------- */
+#include "allocator.hh"
 #include "logger.hh"
 #include "span.hh"
 #include "tamaas.hh"
 #include <memory>
 #include <thrust/copy.h>
 #include <thrust/fill.h>
 #include <utility>
 /* -------------------------------------------------------------------------- */
 
 namespace tamaas {
 
 /// Generic storage class with wrapping capacities
 template <typename T>
 struct Array final {
   /// Default
   Array() = default;
 
   /// Empty array of given size
   Array(UInt size) : Array() { resize(size); }
 
   /// Copy constructor (deep)
   Array(const Array& v) : Array() {
     resize(v.size());
     thrust::copy(v.view_.begin(), v.view_.end(), view_.begin());
   }
 
   /// Move constructor (transfers data ownership)
   Array(Array&& v) noexcept : Array() {
     view_ = std::exchange(v.view_, span<T>{});
     reserved_ = std::exchange(v.reserved_, 0);
     wrapped_ = std::exchange(v.wrapped_, false);
   }
 
   /// Wrap array on data
   Array(T* data, UInt size) noexcept : Array() { wrap(data, size); }
 
   /// Wrap on span
   Array(span<T> view) noexcept : Array() { wrap(view); }
 
   /// Destructor
   ~Array() {
     if (not wrapped_)
       alloc_.deallocate(view_);
   }
 
   /// Copy operator
   Array& operator=(const Array& v) {
     wrapped_ = false;
     resize(v.size());
     thrust::copy(v.view_.begin(), v.view_.end(), view_.begin());
     return *this;
   }
 
   /// Move operator
   Array& operator=(Array&& v) noexcept {
     if (this == &v)
       return *this;
 
     if (not wrapped_)
       alloc_.deallocate(view_);
     view_ = std::exchange(v.view_, span<T>{});
     reserved_ = std::exchange(v.reserved_, 0);
     wrapped_ = std::exchange(v.wrapped_, false);
     return *this;
   }
 
   /// Wrap on view
   Array& operator=(span<T> v) noexcept { wrap(v); }
 
   /// Wrap array
   void wrap(const Array& other) noexcept { wrap(other.view_); }
 
   /// Wrap view
   void wrap(span<T> view) noexcept {
     view_ = view;
     wrapped_ = true;
     reserved_ = 0;
   }
 
   /// Wrap a memory pointer
   void wrap(T* data, UInt size) noexcept { wrap(span<T>{data, size}); }
 
   /// Data pointer access (const)
   const T* data() const { return view_.data(); }
 
   /// Data pointer access (non-const)
   T* data() { return view_.data(); }
 
   /// Resize array
   void resize(UInt new_size, const T& value = T()) {
     if (wrapped_)
       TAMAAS_EXCEPTION("cannot resize wrapped array");
 
     // Erase array
     if (new_size == 0) {
       alloc_.deallocate(view_);
       view_ = span<T>{};
       reserved_ = 0;
       return;
     }
 
     // Do nothing
     if (new_size == size())
       return;
 
     // Allocate new data
     alloc_.deallocate(view_);
 
     view_ = alloc_.allocate(new_size);
     reserved_ = new_size;
 
     if (not wrapped_)
       thrust::fill(view_.begin(), view_.end(), value);
   }
 
   /// Reserve storage space
   void reserve(UInt size) {
     if (reserved_ >= size)
       return;
 
     auto new_view = alloc_.allocate(size);
 
     if (new_view.data() != view_.data()) {
       thrust::copy(view_.begin(), view_.end(), new_view.begin());
 
       alloc_.deallocate(view_);
       view_ = {new_view.data(), view_.size()};
       reserved_ = view_.size();
     } else {
       reserved_ = size;
     }
   }
 
   /// Access operator
   inline T& operator[](UInt i) { return view_[i]; }
 
   /// Access operator (const)
   inline const T& operator[](UInt i) const { return view_[i]; }
 
   /// Get size of array
   inline UInt size() const { return view_.size(); }
 
   span<T> view() const { return view_; }
 
 private:
   span<T> view_;
   typename span<T>::size_type reserved_ = 0;
   bool wrapped_ = false;
   Allocator<T> alloc_;
 };
 
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
 #endif /* ARRAY_HH */
diff --git a/src/core/span.hh b/src/core/span.hh
index a4622de..d5c81d4 100644
--- a/src/core/span.hh
+++ b/src/core/span.hh
@@ -1,71 +1,71 @@
 /*
  *  SPDX-License-Indentifier: AGPL-3.0-or-later
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef SPAN_HH
 #define SPAN_HH
 /* -------------------------------------------------------------------------- */
 #include "tamaas.hh"
 #include <cstddef>
 #include <iterator>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 template <typename T>
 struct span {
   using element_type = T;
   using value_type = std::remove_cv_t<T>;
   using size_type = std::size_t;
   using difference_type = std::ptrdiff_t;
   using pointer = T*;
   using const_pointer = const T*;
   using reference = T&;
   using const_reference = const T&;
   using iterator = T*;
   using reverse_iterator = std::reverse_iterator<iterator>;
 
   constexpr size_type size() const noexcept { return size_; }
   constexpr pointer data() const noexcept { return data_; }
 
   constexpr iterator begin() const noexcept { return data_; }
   constexpr iterator end() const noexcept { return data_ + size_; }
   constexpr iterator begin() noexcept { return data_; }
   constexpr iterator end() noexcept { return data_ + size_; }
 
   reference operator[](size_type idx) {
-    TAMAAS_ASSERT((idx < size_), "index out of span range");
+    TAMAAS_ASSERT(idx < size_, "index out of span range");
     return data_[idx];
   }
 
   const_reference operator[](size_type idx) const {
-    TAMAAS_ASSERT((idx < size_), "index out of span range");
+    TAMAAS_ASSERT(idx < size_, "index out of span range");
     return data_[idx];
   }
 
   pointer data_ = nullptr;
   size_type size_ = 0;
 };
 
 }  // namespace tamaas
 
 /* -------------------------------------------------------------------------- */
 #endif
diff --git a/src/core/tamaas.hh b/src/core/tamaas.hh
index 836c0bf..2af270e 100644
--- a/src/core/tamaas.hh
+++ b/src/core/tamaas.hh
@@ -1,200 +1,186 @@
 /**
  *  @mainpage Tamaas - A high-performance periodic contact library
  *
  *  @section Introduction
  *  Tamaas is a spectral-integral-equation based contact library. It is made
  *  with love to be fast and friendly!
  *
  *  @author Lucas Frérot <lfrerot1@jhu.edu>
  *  @author Guillaume Anciaux <guillaume.anciaux@epfl.ch>
  *  @author Valentine Rey <valentine.rey@univ-nantes.fr>
  *  @author Son Pham-Ba <son.phamba@epfl.ch>
  *  @author Jean-François Molinari <jean-francois.molinari@epfl.ch>
  *
  *  @section License
  *
  *  Copyright (©) 2016-2022 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *  Copyright (©) 2020-2022 Lucas Frérot
  *
  *  This program is free software: you can redistribute it and/or modify
  *  it under the terms of the GNU Affero General Public License as published
  *  by the Free Software Foundation, either version 3 of the License, or
  *  (at your option) any later version.
  *
  *  This program 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 Affero General Public License for more details.
  *
  *  You should have received a copy of the GNU Affero General Public License
  *  along with this program.  If not, see <https://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #ifndef TAMAAS_HH
 #define TAMAAS_HH
 /* -------------------------------------------------------------------------- */
 
 #ifndef TAMAAS_USE_CUDA
 #define TAMAAS_USE_FFTW
 #endif
 
 // Values for fftw backends
 #define TAMAAS_FFTW_BACKEND_OMP 2
 #define TAMAAS_FFTW_BACKEND_THREADS 2
 #define TAMAAS_FFTW_BACKEND_NONE 3
 // Values for thrust backends
 #define TAMAAS_LOOP_BACKEND_OMP 1
 #define TAMAAS_LOOP_BACKEND_TBB 2
 #define TAMAAS_LOOP_BACKEND_CPP 3
 #define TAMAAS_LOOP_BACKEND_CUDA 4
 // Default loop backend is OpenMP
 #ifndef TAMAAS_LOOP_BACKEND
 #define TAMAAS_LOOP_BACKEND TAMAAS_LOOP_BACKEND_OMP
 #endif
 // Default FFTW backend is none
 #ifndef TAMAAS_FFTW_BACKEND
 #define TAMAAS_FFTW_BACKEND TAMAAS_FFTW_BACKEND_NONE
 #endif
 // If the thrust device hasn't been set, set OpenMP
 #ifndef THRUST_DEVICE_SYSTEM
 #define THRUST_DEVICE_SYSTEM THRUST_DEVICE_SYSTEM_OMP
 #endif
 
 /// Convenience macros
 #define TAMAAS_DEBUG_MSG(mesg)                                                 \
   __FILE__ << ':' << __LINE__ << ": " << mesg << '\n'
 #define TAMAAS_EXCEPTION(mesg)                                                 \
   {                                                                            \
     std::stringstream sstr;                                                    \
     sstr << TAMAAS_DEBUG_MSG("FATAL: " << mesg);                               \
     throw ::tamaas::Exception(sstr.str());                                     \
   }
 
 #define SURFACE_FATAL(mesg) TAMAAS_EXCEPTION(mesg)
 
 #if defined(TAMAAS_DEBUG)
 #define TAMAAS_ASSERT(cond, reason)                                            \
   do {                                                                         \
     if (not(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; }
 
 /* -------------------------------------------------------------------------- */
 // Standard includes
 #include <exception>
 #include <iostream>
 #include <memory>
 #include <string>
 #include <type_traits>
 /* -------------------------------------------------------------------------- */
 // Special thrust includes
 #include <thrust/complex.h>
 #include <thrust/random.h>
 
-#ifdef TAMAAS_USE_CUDA
-#include "cuda/unified_allocator.hh"
-#else
-#include "fftw/fftw_allocator.hh"
-#endif
-
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 /* -------------------------------------------------------------------------- */
 
 /// Cuda specific definitions
 #define CUDA_LAMBDA __device__ __host__
 
-#ifdef TAMAAS_USE_CUDA
-template <typename T>
-using Allocator = UnifiedAllocator<T>;
-#else
-template <typename T>
-using Allocator = FFTWAllocator<T>;
-#endif
-
 /// Common types definitions
 // If type macros have not been set, put default values
 #ifndef TAMAAS_REAL_TYPE
 #define TAMAAS_REAL_TYPE double
 #endif
 #ifndef TAMAAS_INT_TYPE
 #define TAMAAS_INT_TYPE int
 #endif
 
 using Real = TAMAAS_REAL_TYPE;           ///< default floating point type
 using Int = TAMAAS_INT_TYPE;             ///< default signed integer type
 using UInt = std::make_unsigned_t<Int>;  ///< default unsigned integer type
 template <typename T>
 using complex = thrust::complex<T>;  ///< template complex wrapper
 using Complex = complex<Real>;       ///< default floating point complex type
 
 /// Defining random toolbox
 using ::thrust::random::normal_distribution;
 using ::thrust::random::uniform_real_distribution;
 using random_engine = ::thrust::random::default_random_engine;
 
 namespace detail {
 template <bool acc, template <typename> class Trait, typename Head,
           typename... Tail>
 struct fold_trait_tail_rec
     : std::integral_constant<bool,
                              fold_trait_tail_rec<acc and Trait<Head>::value,
                                                  Trait, Tail...>::value> {};
 
 template <bool acc, template <typename> class Trait, typename Head>
 struct fold_trait_tail_rec<acc, Trait, Head>
     : std::integral_constant<bool, acc and Trait<Head>::value> {};
 }  // namespace detail
 
 template <template <typename> class Trait, typename... T>
 struct fold_trait : detail::fold_trait_tail_rec<true, Trait, T...> {};
 
 /* -------------------------------------------------------------------------- */
 /// initialize tamaas (0 threads => let OMP_NUM_THREADS decide)
 void initialize(UInt num_threads = 0);
 /// cleanup tamaas
 void finalize();
 /* -------------------------------------------------------------------------- */
 
 /// CUDA-compatible exchange function
 template <class U, class V = U>
 __device__ __host__ U exchange(U& obj, V&& new_value) {
   U old_value = std::move(obj);
   obj = std::forward<V>(new_value);
   return old_value;
 }
 
 /// Generic exception class
 class Exception : public std::exception {
 public:
   /// Constructor
   Exception(std::string mesg) : msg(std::move(mesg)) {}
 
   const char* what() const noexcept override { return msg.c_str(); }
 
   ~Exception() override = default;
 
 private:
   std::string msg;  ///< message of exception
 };
 
 /* -------------------------------------------------------------------------- */
 /// Enumeration of reduction operations
 enum class operation { plus, times, min, max };
 /* -------------------------------------------------------------------------- */
 }  // namespace tamaas
 /* -------------------------------------------------------------------------- */
 
 #endif  // TAMAAS_HH