diff --git a/src/model/integration/accumulator.hh b/src/model/integration/accumulator.hh
index 111dad7..071971c 100644
--- a/src/model/integration/accumulator.hh
+++ b/src/model/integration/accumulator.hh
@@ -1,213 +1,224 @@
 /**
  *  @file
  *  @section LICENSE
  *
  *  Copyright (©) 2016-19 EPFL (École Polytechnique Fédérale de Lausanne),
  *  Laboratory (LSMS - Laboratoire de Simulation en Mécanique des Solides)
  *
  *  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 ACCUMULATOR_HH
 #define ACCUMULATOR_HH
 /* -------------------------------------------------------------------------- */
 #include "grid_hermitian.hh"
 #include "integrator.hh"
 #include "model_type.hh"
 #include "static_types.hh"
+#include <thrust/iterator/zip_iterator.h>
 #include <array>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 /// Accumulation integral manager
 template <model_type type, typename Source,
           typename = std::enable_if_t<is_proxy<Source>::value>>
 class Accumulator {
   using trait = model_type_traits<type>;
   using BufferType = GridHermitian<Real, trait::boundary_dimension>;
 
 public:
   /// Constructor
   Accumulator() {
     // Resizing accumulators (with source)
     for (auto&& acc : accumulator)
       acc.setNbComponents(Source::size);
   }
 
   /// Initialize uniform mesh
   void makeUniformMesh(UInt N, Real domain_size) {
     Real dx = domain_size / (N - 1);
     Real x = 0;
 
     node_positions.resize(N);
     for (UInt i = 0; i < N; ++i, x += dx) {
       node_positions[i] = x;
     }
   }
 
   const std::vector<Real>& nodePositions() const { return node_positions; }
 
   /// Direction flag
   enum class direction { forward, backward };
 
   /// Forward/Backward iterator for integration passes
   template <direction dir>
   struct iterator {
     using integ = Integrator<1>;
     /// Constructor
     iterator(Accumulator& acc, Int k) : acc(acc), k(k) {
       l = (dir != direction::forward) ? k + 1 : k;
     }
     /// Copy constructor
     iterator(const iterator& o)
         : acc(o.acc), k(o.k), l(o.l), r(o.r), xc(o.xc) {}
 
     static constexpr bool upper = dir == direction::backward;
     /// Compare
     bool operator!=(const iterator& o) const { return k != o.k; }
     /// Increment
     iterator& operator++() {
       // update current layer and element. skip last operation
       if (not setup()) {
         next();
         return *this;
       }
 
       Loop::loop(
           [&](auto qv, auto wk0, auto wk1, auto acc_g0, auto acc_g1) {
             Real q = qv.l2norm();
             acc_g0 += wk0 * integ::template G0<upper, 0>(q, r, xc);
             acc_g0 += wk1 * integ::template G0<upper, 1>(q, r, xc);
             acc_g1 += wk0 * integ::template G1<upper, 0>(q, r, xc);
             acc_g1 += wk1 * integ::template G1<upper, 1>(q, r, xc);
           },
           range<VectorProxy<const Real, trait::dimension - 1>>(
               acc.waveVectors())
               .headless(),
           range<Source>(acc.nodeValues()[k]).headless(),
           range<Source>(acc.nodeValues()[k + 1]).headless(),
           range<Source>(acc.accumulator.front()).headless(),
           range<Source>(acc.accumulator.back()).headless());
 
       // accumulating the fundamental mode
       Source wk0(&this->acc.nodeValues()[this->k](0)),
           wk1(&this->acc.nodeValues()[this->k + 1](0));
 
       Source acc_tensor(&acc.accumulator.front()(0));
       acc_tensor += wk0 * integ::template F<0>(this->r);
       acc_tensor += wk1 * integ::template F<1>(this->r);
 
       next();
       return *this;
     }
-    /// Dereference iterator
+    /// Dereference iterator (TODO uniformize tuple types)
     std::tuple<Int, Real, BufferType&, BufferType&> operator*() {
       return {l, acc.node_positions[l], acc.accumulator.front(),
               acc.accumulator.back()};
     }
 
   protected:
     /// Update index layer and element info
     bool setup() {
       if (k + 1 >= static_cast<Int>(acc.node_positions.size()) or k < 0)
         return false;
 
       r = 0.5 * (acc.node_positions[k + 1] - acc.node_positions[k]);
       xc = 0.5 * (acc.node_positions[k + 1] + acc.node_positions[k]);
       return true;
     }
 
     /// Set current layer and update element index
     void next() {
       l = layer(k);
       k = nextElement(k);
     }
 
     /// Element index => layer index
     static Int layer(Int element) {
       return element + (dir == direction::forward);
     }
 
     /// Next element index in right direction
     static Int nextElement(Int element) {
       return element + ((dir == direction::forward) ? 1 : -1);
     }
 
   protected:
     Accumulator& acc;  /// accumulator holder
     Int k = 0;         ///< element index
     Int l = 0;         ///< layer index
     Real r = 0;        ///< element radius
     Real xc = 0;       ///< element center
   };
 
   /// Range for convenience
   template <typename Iterator>
   struct acc_range {
     Iterator _begin, _end;
     Iterator begin() const { return _begin; }
     Iterator end() const { return _end; }
   };
 
   /// Prepare forward loop
   acc_range<iterator<direction::forward>>
   forward(std::vector<BufferType>& nvalues,
           const Grid<Real, trait::boundary_dimension>& wvectors) {
     using it = iterator<direction::forward>;
     reset(nvalues, wvectors);
     return acc_range<it>{it(*this, 0), it(*this, node_positions.size())};
   }
 
   /// Prepare backward loop
   acc_range<iterator<direction::backward>>
   backward(std::vector<BufferType>& nvalues,
            const Grid<Real, trait::boundary_dimension>& wvectors) {
     using it = iterator<direction::backward>;
     reset(nvalues, wvectors);
     // don't ask why the range has those values
     // if it works, don't change it
     return acc_range<it>{it(*this, node_positions.size() - 2), it(*this, -2)};
   }
 
   void reset(std::vector<BufferType>& nvalues,
              const Grid<Real, trait::boundary_dimension>& wvectors) {
     node_values = &nvalues;
     wavevectors = &wvectors;
 
     for (auto&& acc : accumulator) {
       acc.resize(wvectors.sizes());
       acc = 0;
     }
   }
 
   std::vector<BufferType>& nodeValues() {
     TAMAAS_ASSERT(node_values, "Node values is invalid");
     return *node_values;
   }
 
   const Grid<Real, trait::boundary_dimension>& waveVectors() {
     TAMAAS_ASSERT(wavevectors, "Wavevectors is invalid");
     return *wavevectors;
   }
 
+  /// Create a range over the elements in the mesh
+  auto elements() {
+    auto range = Loop::range(node_positions.size());
+    auto begin = thrust::make_zip_iterator(thrust::make_tuple(
+        range.begin(), node_positions.begin(), node_positions.begin() + 1));
+    auto end = thrust::make_zip_iterator(thrust::make_tuple(
+        range.end(), node_positions.end() - 1, node_positions.end()));
+    return acc_range<decltype(begin)>{begin, end};
+  }
+
 private:
   std::array<BufferType, 2> accumulator;
   std::vector<Real> node_positions;
   std::vector<BufferType>* node_values;
   const Grid<Real, trait::boundary_dimension>* wavevectors;
 };
 }  // namespace tamaas
 #endif  // ACCUMULATOR_HH