diff --git a/src/model/integration/accumulator.hh b/src/model/integration/accumulator.hh
index 6887310..c2c5823 100644
--- a/src/model/integration/accumulator.hh
+++ b/src/model/integration/accumulator.hh
@@ -1,242 +1,252 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2018 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 ACCUMULATOR_HH
 #define ACCUMULATOR_HH
 /* -------------------------------------------------------------------------- */
 #include "grid_hermitian.hh"
 #include "integrator.hh"
 #include "model_type.hh"
 #include "static_types.hh"
 #include <array>
 #include <vector>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 
 /// Accumulation integral manager
 template <model_type type, typename Source,
           typename = std::enable_if_t<is_proxy_v<Source>>>
 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;
-  }
+  const std::vector<Real>& nodePositions() const { return node_positions; }
 
   /// Direction flag
   enum class direction { forward, backward };
   /// direction <=> upper
   static constexpr bool is_upper(direction dir) {
     return dir == direction::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(k) {
+    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) {}
     /// 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;
+        next();
+        return *this;
       }
       // part of the domain to be integrated
       constexpr bool upper = is_upper(dir);
 
       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);
+          [&](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()),
-	  range<Source>(acc.nodeValues()[k]),
-	  range<Source>(acc.nodeValues()[k + 1]),
-	  range<Source>(acc.accumulator.front()),
-	  range<Source>(acc.accumulator.back()));
+          },
+          range<VectorProxy<const Real, trait::dimension - 1>>(
+              acc.waveVectors()),
+          range<Source>(acc.nodeValues()[k]),
+          range<Source>(acc.nodeValues()[k + 1]),
+          range<Source>(acc.accumulator.front()),
+          range<Source>(acc.accumulator.back()));
       next();
       return *this;
     }
     /// Dereference iterator
     std::tuple<int, Real, BufferType&, BufferType&> operator*() {
       return {l, acc.node_positions[l], acc.accumulator.front(),
-	      acc.accumulator.back()};
+              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;
+        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 = (dir == direction::forward) ? k + 1 : k;
-      k += (dir == direction::forward) ? 1 : -1;
+      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
+    Int k = 0;         ///< element index
+    Int l = 0;         ///< layer index
+    Real r = 0;        ///< element radius
+    Real xc = 0;       ///< element center
   };
 
   /// Forward/Backward iterator for integration of fundamental mode
   template <direction dir>
   struct fundamental_iterator : iterator<dir> {
     using iterator<dir>::iterator;
     using integ = typename iterator<dir>::integ;
 
     fundamental_iterator& operator++() {
       this->next();
-      //TODO fix
+      // TODO fix
 
       Source wk0(&this->acc.nodeValues()[this->k](0)),
           wk1(&this->acc.nodeValues()[this->k + 1](0));
       acc_tensor += wk0 * integ::template F<0>(this->r);
       acc_tensor += wk1 * integ::template F<1>(this->r);
       return *this;
     }
 
     std::tuple<int, Real, Source> operator*() {
       return {this->l, this->acc.node_positions[this->l], Source(acc_tensor)};
     }
 
   private:
     typename Source::stack_type acc_tensor;
   };
 
   /// 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) {
+          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)};
   }
 
   /// Forward loop for fundamental mode
   acc_range<fundamental_iterator<direction::forward>>
   fundamentalForward(std::vector<BufferType>& nvalues) {
     using it = fundamental_iterator<direction::forward>;
     reset(nvalues, *this->wavevectors);
     return acc_range<it>{it(*this, 0), it(*this, node_positions.size())};
   }
 
   /// Backward loop for fundamental mode
   acc_range<fundamental_iterator<direction::backward>>
   fundamentalBackward(std::vector<BufferType>& nvalues) {
     using it = fundamental_iterator<direction::backward>;
     reset(nvalues, *this->wavevectors);
     return acc_range<it>{it(*this, node_positions.size() - 1), it(*this, -1)};
   }
 
   void reset(std::vector<BufferType>& nvalues,
-	     const Grid<Real, trait::boundary_dimension>& wvectors) {
+             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;
   }
 
 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
diff --git a/src/model/integration/integrator.hh b/src/model/integration/integrator.hh
index b363647..e38729c 100644
--- a/src/model/integration/integrator.hh
+++ b/src/model/integration/integrator.hh
@@ -1,72 +1,73 @@
 /**
  * @file
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2018 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 INTEGRATOR_HH
 #define INTEGRATOR_HH
 /* -------------------------------------------------------------------------- */
 #include "element.hh"
 #include <expolit/expolit>
 /* -------------------------------------------------------------------------- */
 namespace tamaas {
 namespace ex = expolit;
 
 template <UInt interpolation_order>
 class Integrator {
   using element = ExponentialElement<interpolation_order>;
   static constexpr std::pair<Real, Real> bounds{-1, 1};
 
 public:
   /// Standard integral of \f$\exp(\pm qy) \phi(y)\f$ over an element of radius
   /// \f$r\f$ and center \f$x_c\f$
   template <bool upper, UInt shape>
   static Real G0(Real q, Real r, Real xc) {
     const auto F = element::template g0<upper, shape>(q * r);
-    return r * std::exp(q * xc) * ex::definite_integral(bounds, F);
+    return r * std::exp(element::sign(upper) * q * xc) *
+           ex::definite_integral(bounds, F);
   }
 
   /// Standard integral of \f$qy\exp(\pm qy) \phi(y)\f$ over an element of
   /// radius \f$r\f$ and center \f$x_c\f$
   template <bool upper, UInt shape>
   static Real G1(Real q, Real r, Real xc) {
     const auto c = q * r;
     const auto integrals =
         std::make_pair(element::template g0<upper, shape>(c),
                        element::template g1<upper, shape>(c));
     return r * std::exp(element::sign(upper) * q * xc) *
            (q * xc * ex::definite_integral(bounds, integrals.first) +
             ex::definite_integral(bounds, integrals.second));
   }
 
   /// Standard integral of \f$\phi(y)\f$ over an element of radius
   /// \f$r\f$ and center \f$x_c\f$. Used for fundamental mode
   template <UInt shape>
   constexpr static Real F(Real r) {
     return r * ex::definite_integral(bounds, element::shapes[shape]);
   }
 };
 }  // namespace tamaas
 
 #endif  // INTEGRATOR_HH
diff --git a/tests/test_integration.cpp b/tests/test_integration.cpp
index 767080d..0499e1e 100644
--- a/tests/test_integration.cpp
+++ b/tests/test_integration.cpp
@@ -1,147 +1,147 @@
 /**
  *
  * @author Lucas Frérot <lucas.frerot@epfl.ch>
  *
  * @section LICENSE
  *
  * Copyright (©)  2019 EPFL  (Ecole Polytechnique  Fédérale de
  * Lausanne)  Laboratory (LSMS  -  Laboratoire de  Simulation  en Mécanique  des
  * Solides)
  *
  * Tamaas is free  software: you can redistribute it and/or  modify it under the
  * terms  of the  GNU Lesser  General Public  License as  published by  the Free
  * Software Foundation, either version 3 of the License, or (at your option) any
  * later version.
  *
  * Tamaas is  distributed in the  hope that it  will be useful, but  WITHOUT ANY
  * WARRANTY; without even the implied warranty of MERCHANTABILITY or FITNESS FOR
  * A  PARTICULAR PURPOSE. See  the GNU  Lesser General  Public License  for more
  * details.
  *
  * You should  have received  a copy  of the GNU  Lesser General  Public License
  * along with Tamaas. If not, see <http://www.gnu.org/licenses/>.
  *
  */
 /* -------------------------------------------------------------------------- */
 #include "fftransform.hh"
 #include "integration/accumulator.hh"
 #include "test.hh"
 #include <cmath>
 #include <expolit/expolit>
 /* -------------------------------------------------------------------------- */
 using namespace tamaas;
 namespace ex = expolit;
 
 class AccumulatorTest : public ::testing::Test {
 protected:
   void SetUp() override {
     nodal_values.resize(n);
 
     for (auto&& grid : nodal_values) {
       grid.setNbComponents(9);
       grid.resize({2, 2});
       grid = 1;
     }
 
     accumulator.makeUniformMesh(nodal_values.size(), size);
     wavevectors = FFTransform<Real, 2>::computeFrequencies<true>({2, 2});
   }
 
   UInt n = 10;
   Real size = 5;
   std::vector<GridHermitian<Real, 2>> nodal_values;
   Grid<Real, 2> wavevectors;
   Accumulator<model_type::volume_2d, MatrixProxy<Complex, 3, 3>> accumulator;
 };
 
 TEST_F(AccumulatorTest, uniformMesh) {
   auto& pos = this->accumulator.nodePositions();
   std::vector<Real> sol(pos.size());
   std::iota(sol.begin(), sol.end(), 0);
   std::transform(sol.begin(), sol.end(), sol.begin(),
                  [&](auto& x) { return size * x / (sol.size() - 1); });
   ASSERT_TRUE(compare(sol, pos, AreFloatEqual())) << "Uniform mesh fail";
 }
 
 // Litteral type
 using Q = ex::Litteral<struct Q_>;
 
 TEST_F(AccumulatorTest, forward) {
   // Setup for layer checking
   std::vector<int> layers_seen, layers_to_see(n);
   std::iota(layers_to_see.begin(), layers_to_see.end(), 0);
 
   // Setup for integral checking
   constexpr Q q;
   constexpr auto z = ex::Polynomial<Real, 1>({0, 1});
   constexpr auto g0 = ex::exp(q * z);
   constexpr auto g1 = q * z * ex::exp(q * z);
 
   for (auto [l, xl_, acc_g0, acc_g1] :
        accumulator.forward(this->nodal_values, this->wavevectors)) {
     layers_seen.push_back(l);
 
     auto testing = acc_g0(1, 0, 0);
     ASSERT_NEAR(testing.real(), acc_g0(0, 1, 0).real(), 1e-10)
         << "Symmetry fail";
     ASSERT_NEAR(testing.imag(), acc_g0(0, 1, 0).imag(), 1e-10)
         << "Symmetry fail";
 
     testing = acc_g0(1, 1, 0);
     auto ref_value = ex::definite_integral(
         std::make_pair(0., xl_),
         ex::substitute<Q::tag>(ex::Constant<Real>({std::sqrt(2)}), g0));
     ASSERT_NEAR(testing.real(), ref_value, 1e-10)
         << "Integration of exp(qy)*ϕ(y) fail";
 
     testing = acc_g1(1, 1, 0);
     ref_value = ex::definite_integral(
         std::make_pair(0., xl_),
         ex::substitute<Q::tag>(ex::Constant<Real>({std::sqrt(2)}), g1));
     ASSERT_NEAR(testing.real(), ref_value, 1e-10)
         << "Integration of qy*exp(qy)*ϕ(y) fail";
   }
 
   ASSERT_TRUE(compare(layers_seen, layers_to_see))
       << "Did not see correct layers";
 }
 
 TEST_F(AccumulatorTest, backward) {
   // Setup for layer checking
   std::vector<int> layers_seen, layers_to_see(n);
   std::iota(layers_to_see.rbegin(), layers_to_see.rend(), 0);
 
   // Setup for integral checking
   constexpr Q q;
   constexpr auto z = ex::Polynomial<Real, 1>({0, 1});
   constexpr auto g0 = ex::exp(q * (z * (-1)));
   constexpr auto g1 = q * z * ex::exp(q * ((-1) * z));
 
   for (auto [l, xl_, acc_g0, acc_g1] :
        accumulator.backward(nodal_values, wavevectors)) {
     layers_seen.push_back(l);
 
     auto testing = acc_g0(1, 0, 0);
     ASSERT_NEAR(testing.real(), acc_g0(0, 1, 0).real(), 1e-10)
         << "Symmetry fail";
     ASSERT_NEAR(testing.imag(), acc_g0(0, 1, 0).imag(), 1e-10)
         << "Symmetry fail";
 
     testing = acc_g0(1, 1, 0);
     auto ref_value = ex::definite_integral(
         std::make_pair(xl_, size),
         ex::substitute<Q::tag>(ex::Constant<Real>({std::sqrt(2)}), g0));
-    ASSERT_NEAR(testing.real(), ref_value, 1e-10)
+    EXPECT_NEAR(testing.real(), ref_value, 1e-10)
         << "Integration of exp(-qy)*ϕ(y) fail";
 
     testing = acc_g1(1, 1, 0);
     ref_value = ex::definite_integral(
         std::make_pair(xl_, size),
         ex::substitute<Q::tag>(ex::Constant<Real>({std::sqrt(2)}), g1));
-    ASSERT_NEAR(testing.real(), ref_value, 1e-10)
+    EXPECT_NEAR(testing.real(), ref_value, 1e-10)
         << "Integration of qy*exp(-qy)*ϕ(y) fail";
   }
 
   ASSERT_TRUE(compare(layers_seen, layers_to_see))
       << "Did not see correct layers";
 }