diff --git a/src/include/psc/current_deposition.hxx b/src/include/psc/current_deposition.hxx
index f46d5a161c..4447bbb71c 100644
--- a/src/include/psc/current_deposition.hxx
+++ b/src/include/psc/current_deposition.hxx
@@ -10,93 +10,79 @@ class CurrentDeposition1vb
 public:
   using Curr = C;
   using real_t = typename Curr::real_t;
-  using real3_t = Vec3<real_t>;
+  using Real3 = Vec3<real_t>;
 
-  GT_INLINE CurrentDeposition1vb(real3_t fnqs) : fnqs_{fnqs} {}
+  GT_INLINE CurrentDeposition1vb(Real3 fnqs) : fnqs_{fnqs} {}
 
-  GT_INLINE void operator()(Curr& curr, const int i[3], real_t qni_wni,
-                            const real_t dx[3], const real_t xa[3],
+  GT_INLINE void operator()(Curr& curr, const Int3& i, real_t qni_wni,
+                            const Real3& dx, const Real3& xa,
                             dim_xyz tag_dim) const
   {
-    real_t h = (1.f / real_t(12.f)) * dx[0] * dx[1] * dx[2];
-    real3_t fnq = qni_wni * fnqs_;
-
-    curr.add(0, i[0], i[1], i[2],
-             fnq[0] * (dx[0] * (1.f - xa[1]) * (1.f - xa[2]) + h));
-    curr.add(0, i[0], i[1] + 1, i[2],
-             fnq[0] * (dx[0] * (xa[1]) * (1.f - xa[2]) - h));
-    curr.add(0, i[0], i[1], i[2] + 1,
-             fnq[0] * (dx[0] * (1.f - xa[1]) * (xa[2]) - h));
-    curr.add(0, i[0], i[1] + 1, i[2] + 1,
-             fnq[0] * (dx[0] * (xa[1]) * (xa[2]) + h));
-
-    curr.add(1, i[0], i[1], i[2],
-             fnq[1] * (dx[1] * (1.f - xa[0]) * (1.f - xa[2]) + h));
-    curr.add(1, i[0] + 1, i[1], i[2],
-             fnq[1] * (dx[1] * (xa[0]) * (1.f - xa[2]) - h));
-    curr.add(1, i[0], i[1], i[2] + 1,
-             fnq[1] * (dx[1] * (1.f - xa[0]) * (xa[2]) - h));
-    curr.add(1, i[0] + 1, i[1], i[2] + 1,
-             fnq[1] * (dx[1] * (xa[0]) * (xa[2]) + h));
-
-    curr.add(2, i[0], i[1], i[2],
-             fnq[2] * (dx[2] * (1.f - xa[0]) * (1.f - xa[1]) + h));
-    curr.add(2, i[0] + 1, i[1], i[2],
-             fnq[2] * (dx[2] * (xa[0]) * (1.f - xa[1]) - h));
-    curr.add(2, i[0], i[1] + 1, i[2],
-             fnq[2] * (dx[2] * (1.f - xa[0]) * (xa[1]) - h));
-    curr.add(2, i[0] + 1, i[1] + 1, i[2],
-             fnq[2] * (dx[2] * (xa[0]) * (xa[1]) + h));
+    real_t h = (1.f / real_t(12.f)) * dx.prod();
+    Real3 fnq = qni_wni * fnqs_;
+
+    // clang-format off
+    curr.add(0, i[0],     i[1],     i[2],     fnq[0] * (dx[0] * (1.f - xa[1]) * (1.f - xa[2]) + h));
+    curr.add(0, i[0],     i[1] + 1, i[2],     fnq[0] * (dx[0] * (      xa[1]) * (1.f - xa[2]) - h));
+    curr.add(0, i[0],     i[1],     i[2] + 1, fnq[0] * (dx[0] * (1.f - xa[1]) * (      xa[2]) - h));
+    curr.add(0, i[0],     i[1] + 1, i[2] + 1, fnq[0] * (dx[0] * (      xa[1]) * (      xa[2]) + h));
+
+    curr.add(1, i[0],     i[1],     i[2],     fnq[1] * (dx[1] * (1.f - xa[2]) * (1.f - xa[0]) + h));
+    curr.add(1, i[0],     i[1],     i[2] + 1, fnq[1] * (dx[1] * (      xa[2]) * (1.f - xa[0]) - h));
+    curr.add(1, i[0] + 1, i[1],     i[2],     fnq[1] * (dx[1] * (1.f - xa[2]) * (      xa[0]) - h));
+    curr.add(1, i[0] + 1, i[1],     i[2] + 1, fnq[1] * (dx[1] * (      xa[2]) * (      xa[0]) + h));
+
+    curr.add(2, i[0],     i[1],     i[2],     fnq[2] * (dx[2] * (1.f - xa[0]) * (1.f - xa[1]) + h));
+    curr.add(2, i[0] + 1, i[1],     i[2],     fnq[2] * (dx[2] * (      xa[0]) * (1.f - xa[1]) - h));
+    curr.add(2, i[0],     i[1] + 1, i[2],     fnq[2] * (dx[2] * (1.f - xa[0]) * (      xa[1]) - h));
+    curr.add(2, i[0] + 1, i[1] + 1, i[2],     fnq[2] * (dx[2] * (      xa[0]) * (      xa[1]) + h));
+    // clang-format on
   };
 
-  GT_INLINE void operator()(Curr& curr, const int i[3], real_t qni_wni,
-                            const real_t dx[3], const real_t xa[3],
+  GT_INLINE void operator()(Curr& curr, const Int3& i, real_t qni_wni,
+                            const Real3& dx, const Real3& xa,
                             dim_xz tag_dim) const
   {
-    real_t h = (1.f / real_t(12.f)) * dx[0] * dx[1] * dx[2];
-    real3_t fnq = qni_wni * fnqs_;
-
-    curr.add(0, i[0], i[1], i[2], fnq[0] * (dx[0] * (1.f - xa[2])));
-    curr.add(0, i[0], i[1], i[2] + 1, fnq[0] * (dx[0] * (xa[2])));
-
-    curr.add(1, i[0], i[1], i[2],
-             fnq[1] * (dx[1] * (1.f - xa[0]) * (1.f - xa[2]) + h));
-    curr.add(1, i[0] + 1, i[1], i[2],
-             fnq[1] * (dx[1] * (xa[0]) * (1.f - xa[2]) - h));
-    curr.add(1, i[0], i[1], i[2] + 1,
-             fnq[1] * (dx[1] * (1.f - xa[0]) * (xa[2]) - h));
-    curr.add(1, i[0] + 1, i[1], i[2] + 1,
-             fnq[1] * (dx[1] * (xa[0]) * (xa[2]) + h));
-
-    curr.add(2, i[0], i[1], i[2], fnq[2] * (dx[2] * (1.f - xa[0])));
-    curr.add(2, i[0] + 1, i[1], i[2], fnq[2] * (dx[2] * (xa[0])));
+    real_t h = (1.f / real_t(12.f)) * dx.prod();
+    Real3 fnq = qni_wni * fnqs_;
+
+    // clang-format off
+    curr.add(0, i[0],     i[1],     i[2],     fnq[0] * (dx[0] *                 (1.f - xa[2])    ));
+    curr.add(0, i[0],     i[1],     i[2] + 1, fnq[0] * (dx[0] *                 (      xa[2])    ));
+
+    curr.add(1, i[0],     i[1],     i[2],     fnq[1] * (dx[1] * (1.f - xa[2]) * (1.f - xa[0]) + h));
+    curr.add(1, i[0],     i[1],     i[2] + 1, fnq[1] * (dx[1] * (      xa[2]) * (1.f - xa[0]) - h));
+    curr.add(1, i[0] + 1, i[1],     i[2],     fnq[1] * (dx[1] * (1.f - xa[2]) * (      xa[0]) - h));
+    curr.add(1, i[0] + 1, i[1],     i[2] + 1, fnq[1] * (dx[1] * (      xa[2]) * (      xa[0]) + h));
+
+    curr.add(2, i[0],     i[1],     i[2],     fnq[2] * (dx[2] * (1.f - xa[0])                    ));
+    curr.add(2, i[0] + 1, i[1],     i[2],     fnq[2] * (dx[2] * (      xa[0])                    ));
+    // clang-format on
   }
 
-  GT_INLINE void operator()(Curr& curr, const int i[3], real_t qni_wni,
-                            const real_t dx[3], const real_t xa[3],
+  GT_INLINE void operator()(Curr& curr, const Int3& i, real_t qni_wni,
+                            const Real3& dx, const Real3& xa,
                             dim_yz tag_dim) const
   {
-    real_t h = (1.f / real_t(12.f)) * dx[0] * dx[1] * dx[2];
-
-    real3_t fnq = qni_wni * fnqs_;
-    curr.add(0, i[0], i[1], i[2],
-             fnq[0] * (dx[0] * (1.f - xa[1]) * (1.f - xa[2]) + h));
-    curr.add(0, i[0], i[1] + 1, i[2],
-             fnq[0] * (dx[0] * (xa[1]) * (1.f - xa[2]) - h));
-    curr.add(0, i[0], i[1], i[2] + 1,
-             fnq[0] * (dx[0] * (1.f - xa[1]) * (xa[2]) - h));
-    curr.add(0, i[0], i[1] + 1, i[2] + 1,
-             fnq[0] * (dx[0] * (xa[1]) * (xa[2]) + h));
-
-    curr.add(1, i[0], i[1], i[2], fnq[1] * (dx[1] * (1.f - xa[2])));
-    curr.add(1, i[0], i[1], i[2] + 1, fnq[1] * (dx[1] * (xa[2])));
-
-    curr.add(2, i[0], i[1], i[2], fnq[2] * (dx[2] * (1.f - xa[1])));
-    curr.add(2, i[0], i[1] + 1, i[2], fnq[2] * (dx[2] * (xa[1])));
+    real_t h = (1.f / real_t(12.f)) * dx.prod();
+    Real3 fnq = qni_wni * fnqs_;
+
+    // clang-format off
+    curr.add(0, i[0],     i[1],     i[2],     fnq[0] * (dx[0] * (1.f - xa[1]) * (1.f - xa[2]) + h));
+    curr.add(0, i[0],     i[1] + 1, i[2],     fnq[0] * (dx[0] * (      xa[1]) * (1.f - xa[2]) - h));
+    curr.add(0, i[0],     i[1],     i[2] + 1, fnq[0] * (dx[0] * (1.f - xa[1]) * (      xa[2]) - h));
+    curr.add(0, i[0],     i[1] + 1, i[2] + 1, fnq[0] * (dx[0] * (      xa[1]) * (      xa[2]) + h));
+
+    curr.add(1, i[0],     i[1],     i[2],     fnq[1] * (dx[1] * (1.f - xa[2])                    ));
+    curr.add(1, i[0],     i[1],     i[2] + 1, fnq[1] * (dx[1] * (      xa[2])                    ));
+
+    curr.add(2, i[0],     i[1],     i[2],     fnq[2] * (dx[2] *                 (1.f - xa[1])    ));
+    curr.add(2, i[0],     i[1] + 1, i[2],     fnq[2] * (dx[2] *                 (      xa[1])    ));
+    // clang-format on
   }
 
 private:
-  real3_t fnqs_;
+  Real3 fnqs_;
 };
 
 }; // namespace psc
diff --git a/src/include/psc/deposit.hxx b/src/include/psc/deposit.hxx
index af0a89f3f8..f00d3a9235 100644
--- a/src/include/psc/deposit.hxx
+++ b/src/include/psc/deposit.hxx
@@ -26,13 +26,13 @@ class Deposit1st
 {
 public:
   using real_t = R;
-  using real3_t = gt::sarray<R, 3>;
+  using Real3 = Vec3<real_t>;
   using dim_t = D;
 
 private:
   template <typename F>
-  void operator()(F& flds, const gt::sarray<int, 3>& l, const real3_t& h,
-                  real_t val, dim_yz tag)
+  void operator()(F& flds, const Int3& l, const Real3& h, real_t val,
+                  dim_yz tag)
   {
     flds(0, l[1] + 0, l[2] + 0) += val * (1.f - h[1]) * (1.f - h[2]);
     flds(0, l[1] + 1, l[2] + 0) += val * h[1] * (1.f - h[2]);
@@ -41,8 +41,8 @@ private:
   }
 
   template <typename F>
-  void operator()(F& flds, const gt::sarray<int, 3>& l, const real3_t& h,
-                  real_t val, dim_xyz tag)
+  void operator()(F& flds, const Int3& l, const Real3& h, real_t val,
+                  dim_xyz tag)
   {
     // clang-format off
     flds(l[0] + 0, l[1] + 0, l[2] + 0) += val * (1.f - h[0]) * (1.f - h[1]) * (1.f - h[2]);
@@ -58,8 +58,7 @@ private:
 
 public:
   template <typename F>
-  void operator()(F& flds, const gt::sarray<int, 3>& l, const real3_t& h,
-                  real_t val)
+  void operator()(F& flds, const Int3& l, const Real3& h, real_t val)
   {
     (*this)(flds, l, h, val, dim_t{});
   }
@@ -76,13 +75,13 @@ class Deposit2nd
 {
 public:
   using real_t = R;
-  using real3_t = gt::sarray<R, 3>;
+  using Real3 = Vec3<R>;
   using dim_t = D;
 
 private:
   template <typename F>
-  void operator()(F& flds, const gt::sarray<int, 3>& l, const real3_t& h,
-                  real_t val, dim_yz tag)
+  void operator()(F& flds, const Int3& l, const Real3& h, real_t val,
+                  dim_yz tag)
   {
     real_t gmy = .5f * (.5f + h[1]) * (.5f + h[1]);
     real_t gmz = .5f * (.5f + h[2]) * (.5f + h[2]);
@@ -105,8 +104,8 @@ private:
   }
 
   template <typename F>
-  void operator()(F& flds, const gt::sarray<int, 3>& l, const real3_t& h,
-                  real_t val, dim_xyz tag)
+  void operator()(F& flds, const Int3& l, const Real3& h, real_t val,
+                  dim_xyz tag)
   {
     real_t gmx = .5f * (.5f + h[0]) * (.5f + h[0]);
     real_t gmy = .5f * (.5f + h[1]) * (.5f + h[1]);
@@ -151,8 +150,7 @@ private:
 
 public:
   template <typename F>
-  void operator()(F& flds, const gt::sarray<int, 3>& l, const real3_t& h,
-                  real_t val)
+  void operator()(F& flds, const Int3& l, const Real3& h, real_t val)
   {
     (*this)(flds, l, h, val, dim_t{});
   }
@@ -178,12 +176,10 @@ public:
   static std::string suffix() { return "_1st_nc"; }
 
   template <typename F>
-  void operator()(F& flds, const gt::sarray<int, 3>& ib,
-                  const gt::sarray<T, 3>& x, T val)
-
+  void operator()(F& flds, const Int3& ib, const Vec3<T>& x, T val)
   {
-    gt::sarray<int, 3> l;
-    gt::sarray<T, 3> h;
+    Int3 l;
+    Vec3<T> h;
     for (int d = 0; d < 3; d++) {
       l[d] = fint(x[d]);
       h[d] = x[d] - l[d];
@@ -201,12 +197,10 @@ public:
   static std::string suffix() { return "_1st_cc"; }
 
   template <typename F>
-  void operator()(F& flds, const gt::sarray<int, 3>& ib,
-                  const gt::sarray<T, 3>& x, T val)
-
+  void operator()(F& flds, const Int3& ib, const Vec3<T>& x, T val)
   {
-    gt::sarray<int, 3> l;
-    gt::sarray<T, 3> h;
+    Int3 l;
+    Vec3<T> h;
     for (int d = 0; d < 3; d++) {
       l[d] = fint(x[d] - .5f);
       h[d] = x[d] - .5f - l[d];
@@ -228,12 +222,10 @@ public:
   static std::string suffix() { return "_2nd_nc"; }
 
   template <typename F>
-  void operator()(F& flds, const gt::sarray<int, 3>& ib,
-                  const gt::sarray<T, 3>& x, T val)
-
+  void operator()(F& flds, const Int3& ib, const Vec3<T>& x, T val)
   {
-    gt::sarray<int, 3> l;
-    gt::sarray<T, 3> h;
+    Int3 l;
+    Vec3<T> h;
     for (int d = 0; d < 3; d++) {
       l[d] = nint(x[d]);
       h[d] = l[d] - x[d]; // negated!
@@ -245,14 +237,14 @@ public:
 };
 
 template <typename D, typename F, typename T>
-void nc(F& flds, const gt::sarray<int, 3>& ib, const gt::sarray<T, 3>& x, T val)
+void nc(F& flds, const Int3& ib, const Vec3<T>& x, T val)
 {
   Deposit1stNc<T, D> deposit;
   deposit(flds, ib, x, val);
 }
 
 template <typename D, typename F, typename T>
-void cc(F& flds, const gt::sarray<int, 3>& ib, const gt::sarray<T, 3>& x, T val)
+void cc(F& flds, const Int3& ib, const Vec3<T>& x, T val)
 {
   Deposit1stCc<T, D> deposit;
   deposit(flds, ib, x, val);
@@ -276,27 +268,25 @@ class Deposit
 public:
   using real_t = R;
   using dim_t = D;
-  using real3_t = gt::sarray<real_t, 3>;
+  using Real3 = Vec3<real_t>;
   using DepositNorm = DEPOSIT_NORM<real_t, dim_t>;
   static const centering::Centering CENTERING = DepositNorm::CENTERING;
 
   static std::string suffix() { return DepositNorm::suffix(); }
 
-  Deposit(const real3_t& dx, real_t fnqs) : dxi_{real_t(1.) / dx}, fnqs_{fnqs}
-  {}
+  Deposit(const Real3& dx, real_t fnqs) : dxi_{real_t(1.) / dx}, fnqs_{fnqs} {}
 
   template <typename F>
-  void operator()(const F& flds, const gt::shape_type<3>& ib, const real3_t& xi,
-                  real_t val)
+  void operator()(const F& flds, const Int3& ib, const Real3& xi, real_t val)
   {
-    real3_t x = xi * dxi_;
+    Real3 x = xi * dxi_;
     real_t value = fnqs_ * val;
 
     DepositNorm deposit;
     deposit(flds, ib, x, value);
   }
 
-  real3_t dxi_;
+  Real3 dxi_;
   real_t fnqs_;
 };
 
diff --git a/src/include/psc/gtensor.h b/src/include/psc/gtensor.h
index 98ec8ebb61..7ca29bdb4c 100644
--- a/src/include/psc/gtensor.h
+++ b/src/include/psc/gtensor.h
@@ -8,92 +8,3 @@
 #endif
 
 #include <gtensor/gtensor.h>
-
-namespace
-{
-
-template <typename T, gt::size_type N>
-GT_INLINE gt::sarray<T, N> operator*(T a, const gt::sarray<T, N>& x)
-{
-  gt::sarray<T, N> res;
-  for (gt::size_type i = 0; i < N; i++) {
-    res[i] = a * x[i];
-  }
-  return res;
-}
-
-template <typename T, gt::size_type N>
-GT_INLINE gt::sarray<T, N> operator/(const gt::sarray<T, N>& x, T a)
-{
-  gt::sarray<T, N> res;
-  for (gt::size_type i = 0; i < N; i++) {
-    res[i] = x[i] / a;
-  }
-  return res;
-}
-
-template <typename T, gt::size_type N>
-GT_INLINE gt::sarray<T, N> operator/(T a, const gt::sarray<T, N>& x)
-{
-  gt::sarray<T, N> res;
-  for (gt::size_type i = 0; i < N; i++) {
-    res[i] = a / x[i];
-  }
-  return res;
-}
-
-template <typename T, gt::size_type N>
-GT_INLINE gt::sarray<T, N> operator-(const gt::sarray<T, N>& x)
-{
-  gt::sarray<T, N> res;
-  for (gt::size_type i = 0; i < N; i++) {
-    res[i] = -x[i];
-  }
-  return res;
-}
-
-template <typename T, gt::size_type N>
-GT_INLINE gt::sarray<T, N> operator+(const gt::sarray<T, N>& x,
-                                     const gt::sarray<T, N>& y)
-{
-  gt::sarray<T, N> res;
-  for (gt::size_type i = 0; i < N; i++) {
-    res[i] = x[i] + y[i];
-  }
-  return res;
-}
-
-template <typename T, gt::size_type N>
-GT_INLINE gt::sarray<T, N> operator-(const gt::sarray<T, N>& x,
-                                     const gt::sarray<T, N>& y)
-{
-  gt::sarray<T, N> res;
-  for (gt::size_type i = 0; i < N; i++) {
-    res[i] = x[i] - y[i];
-  }
-  return res;
-}
-
-template <typename T, gt::size_type N>
-GT_INLINE gt::sarray<T, N> operator*(const gt::sarray<T, N>& x,
-                                     const gt::sarray<T, N>& y)
-{
-  gt::sarray<T, N> res;
-  for (gt::size_type i = 0; i < N; i++) {
-    res[i] = x[i] * y[i];
-  }
-  return res;
-}
-
-template <typename T, gt::size_type N>
-GT_INLINE gt::sarray<T, N> operator/(const gt::sarray<T, N>& x,
-                                     const gt::sarray<T, N>& y)
-{
-  gt::sarray<T, N> res;
-  for (gt::size_type i = 0; i < N; i++) {
-    res[i] = x[i] / y[i];
-  }
-  return res;
-}
-
-} // namespace
diff --git a/src/include/psc/moment.hxx b/src/include/psc/moment.hxx
index 39672f7ed5..08a108714d 100644
--- a/src/include/psc/moment.hxx
+++ b/src/include/psc/moment.hxx
@@ -18,11 +18,11 @@ struct DepositParticlesContext
   using storage_type = S;
   using dim_t = D;
   using real_t = typename storage_type::value_type;
-  using real3_t = gt::sarray<real_t, 3>;
+  using Real3 = Vec3<real_t>;
   using Deposit = DepositCode<real_t, dim_t>;
 
   DepositParticlesContext(storage_type& mflds_gt, const Int3& ib,
-                          const real3_t& dx, real_t fnqs)
+                          const Real3& dx, real_t fnqs)
     : mflds_gt{mflds_gt}, ib{ib}, deposit{dx, fnqs}
   {}
 
@@ -35,7 +35,7 @@ struct DepositParticlesContext
   Int3 ib;
   Deposit deposit;
   int p;
-  real3_t x;
+  Real3 x;
 };
 
 template <typename S, typename D,
@@ -46,11 +46,11 @@ public:
   using storage_type = S;
   using dim_t = D;
   using real_t = typename storage_type::value_type;
-  using real3_t = gt::sarray<real_t, 3>;
+  using Real3 = Vec3<real_t>;
   using DepositCtx = DepositParticlesContext<storage_type, dim_t, DepositCode>;
 
   template <typename Mparticles, typename F>
-  void operator()(const real3_t& dx, real_t fnqs, storage_type& mflds_gt,
+  void operator()(const Real3& dx, real_t fnqs, storage_type& mflds_gt,
                   const Int3& ib, const Mparticles& mprts, F&& func)
   {
     auto accessor = mprts.accessor();
@@ -70,10 +70,10 @@ template <template <typename, typename> class DepositCode, typename D,
 void deposit(MF& mflds_gt, const Int3& ib, const MP& mprts, F&& func)
 {
   using real_t = typename MF::value_type;
-  using real3_t = gt::sarray<real_t, 3>;
+  using Real3 = Vec3<real_t>;
 
   const auto& domain = mprts.grid().domain;
-  real3_t dx = {domain.dx[0], domain.dx[1], domain.dx[2]};
+  Real3 dx = {domain.dx[0], domain.dx[1], domain.dx[2]};
   real_t fnqs = mprts.grid().norm.fnqs;
   DepositParticles<MF, D, DepositCode> deposit;
   deposit(dx, fnqs, mflds_gt, ib, mprts, std::forward<F>(func));
diff --git a/src/libpsc/psc_push_fields/marder_impl.hxx b/src/libpsc/psc_push_fields/marder_impl.hxx
index 436c0c7bb5..7acb5696fc 100644
--- a/src/libpsc/psc_push_fields/marder_impl.hxx
+++ b/src/libpsc/psc_push_fields/marder_impl.hxx
@@ -61,14 +61,14 @@ inline void correct(const Grid_t& grid, E1& efield, const Int3& efield_ib,
                     typename E1::value_type diffusion)
 {
   using real_t = gt::expr_value_type<E1>;
-  using real3_t = gt::sarray<real_t, 3>;
+  using Real3 = Vec3<real_t>;
 
   assert(efield_ib == -grid.ibn);
   assert(mf_ib == -grid.ibn);
 
-  real3_t fac = {.5f * grid.dt * diffusion / grid.domain.dx[0],
-                 .5f * grid.dt * diffusion / grid.domain.dx[1],
-                 .5f * grid.dt * diffusion / grid.domain.dx[2]};
+  Real3 fac = {.5f * grid.dt * diffusion / grid.domain.dx[0],
+               .5f * grid.dt * diffusion / grid.domain.dx[1],
+               .5f * grid.dt * diffusion / grid.domain.dx[2]};
 
   for (int p = 0; p < grid.n_patches(); p++) {
     Int3 lx, rx, ly, ry, lz, rz;
diff --git a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx
index 393ed033bf..c6af987484 100644
--- a/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx
+++ b/src/libpsc/psc_push_particles/inc_curr_1vb_split.cxx
@@ -1,6 +1,9 @@
 
 #pragma once
 
+// J. Villasenor and O. Buneman, "Rigorous charge conservation for local
+// electromagnetic field solvers", Computer Physics Communications 69 (1992) 306
+
 // ======================================================================
 
 template <typename Order, typename Dim, typename _fields_t>
@@ -15,40 +18,28 @@ struct Current1vbSplit
 
   Current1vbSplit(const Grid_t& grid)
     : dt_(grid.dt),
-      dxi_{Real3{1., 1., 1.} / Real3(grid.domain.dx)},
+      dxi_{grid.domain.dx.inv()},
       deposition_(real_t(grid.norm.fnqs / grid.dt) * Real3(grid.domain.dx))
   {}
 
-  void calc_j2_one_cell(fields_t curr_cache, real_t qni_wni, const real_t xm[3],
-                        const real_t xp[3])
+  void calc_j2_one_cell(fields_t curr_cache, real_t qni_wni, const Real3& xm,
+                        const Real3& xp)
   {
-    real_t dx[3] = {xp[0] - xm[0], xp[1] - xm[1], xp[2] - xm[2]};
-    real_t xa[3] = {.5f * (xm[0] + xp[0]), .5f * (xm[1] + xp[1]),
-                    .5f * (xm[2] + xp[2])};
-
-    int i[3];
-    for (int d = 0; d < 3; d++) {
-      i[d] = fint(xa[d]);
-      xa[d] -= i[d];
-    }
-
+    Real3 dx = xp - xm;
+    Real3 xa = real_t(0.5) * (xp + xm);
+    Int3 i = xa.fint();
+    xa -= Real3(i);
     deposition_(curr_cache, i, qni_wni, dx, xa, Dim{});
   }
 
-  static void calc_j2_split_along_dim(int dim, int im, real_t x1[3],
-                                      const real_t xm[3], const real_t xp[3])
+  static Real3 calc_split_x1(int dim, int im, int ip, const Real3& xm,
+                             const Real3& xp)
   {
-    real_t bnd = 0.f; // quell warning
-    int ip = fint(xp[dim]);
-    if (ip == im + 1) { // crossed boundary to right
-      bnd = ip;
-    } else if (ip == im - 1) { // crosses boundary to left
-      bnd = im;
-    } else {
-      assert(0);
-    }
+    // boundary is at the lower edge of the cell with the higher index
+    real_t bnd = std::max(im, ip);
     real_t frac = (bnd - xm[dim]) / (xp[dim] - xm[dim]);
-    // FIXME, set d == dim value to exact boundary?
+
+    Real3 x1; // = xm + frac * (xp - xm), but want exact x1[d]
     for (int d = 0; d < 3; d++) {
       if (d == dim) {
         x1[d] = bnd;
@@ -56,58 +47,61 @@ struct Current1vbSplit
         x1[d] = xm[d] + frac * (xp[d] - xm[d]);
       }
     }
+    return x1;
   }
 
-  void calc_j2_split_dim_x(fields_t curr_cache, real_t qni_wni,
-                           const real_t* xm, const real_t* xp)
+  // TODO c++17 combine these calc_j2_split_dim_* with if constexpr
+
+  void calc_j2_split_dim_x(fields_t curr_cache, real_t qni_wni, const Real3& xm,
+                           const Real3& xp)
   {
-    const int dim = 0;
+    constexpr int dim = 0;
     int im = fint(xm[dim]);
-    if (!Dim::InvarX::value && (fint(xp[dim]) != im)) {
-      real_t x1[3];
-      calc_j2_split_along_dim(dim, im, x1, xm, xp);
+    int ip = fint(xp[dim]);
+    if (Dim::is_invar(dim) || ip == im) {
+      calc_j2_one_cell(curr_cache, qni_wni, xm, xp);
+    } else {
+      Real3 x1 = calc_split_x1(dim, im, ip, xm, xp);
       calc_j2_one_cell(curr_cache, qni_wni, xm, x1);
       calc_j2_one_cell(curr_cache, qni_wni, x1, xp);
-    } else {
-      calc_j2_one_cell(curr_cache, qni_wni, xm, xp);
     }
   }
 
-  void calc_j2_split_dim_y(fields_t curr_cache, real_t qni_wni,
-                           const real_t* xm, const real_t* xp)
+  void calc_j2_split_dim_y(fields_t curr_cache, real_t qni_wni, const Real3& xm,
+                           const Real3& xp)
   {
-    const int dim = 1;
+    constexpr int dim = 1;
     int im = fint(xm[dim]);
-    if (!Dim::InvarY::value && (fint(xp[dim]) != im)) {
-      real_t x1[3];
-      calc_j2_split_along_dim(dim, im, x1, xm, xp);
+    int ip = fint(xp[dim]);
+    if (Dim::is_invar(dim) || ip == im) {
+      calc_j2_split_dim_x(curr_cache, qni_wni, xm, xp);
+    } else {
+      Real3 x1 = calc_split_x1(dim, im, ip, xm, xp);
       calc_j2_split_dim_x(curr_cache, qni_wni, xm, x1);
       calc_j2_split_dim_x(curr_cache, qni_wni, x1, xp);
-    } else {
-      calc_j2_split_dim_x(curr_cache, qni_wni, xm, xp);
     }
   }
 
-  void calc_j2_split_dim_z(fields_t curr_cache, real_t qni_wni,
-                           const real_t* xm, const real_t* xp)
+  void calc_j2_split_dim_z(fields_t curr_cache, real_t qni_wni, const Real3& xm,
+                           const Real3& xp)
   {
-    const int dim = 2;
+    constexpr int dim = 2;
     int im = fint(xm[dim]);
-    if (!Dim::InvarZ::value && (fint(xp[dim]) != im)) {
-      real_t x1[3];
-      calc_j2_split_along_dim(dim, im, x1, xm, xp);
+    int ip = fint(xp[dim]);
+    if (Dim::is_invar(dim) || ip == im) {
+      calc_j2_split_dim_y(curr_cache, qni_wni, xm, xp);
+    } else {
+      Real3 x1 = calc_split_x1(dim, im, ip, xm, xp);
       calc_j2_split_dim_y(curr_cache, qni_wni, xm, x1);
       calc_j2_split_dim_y(curr_cache, qni_wni, x1, xp);
-    } else {
-      calc_j2_split_dim_y(curr_cache, qni_wni, xm, xp);
     }
   }
 
   // ----------------------------------------------------------------------
   // calc_j
 
-  void calc_j(fields_t curr_cache, real_t* xm, real_t* xp, int* lf, int* lg,
-              real_t qni_wni, real_t* vxi)
+  void calc_j(fields_t curr_cache, Real3& xm, Real3& xp, const Int3& lf,
+              const Int3& lg, real_t qni_wni, Real3& vxi)
   {
     // this way, we guarantee that the average position in invar directions
     // will remain in the 0th cell
diff --git a/src/libpsc/psc_push_particles/inc_curr_zigzag.cxx b/src/libpsc/psc_push_particles/inc_curr_zigzag.cxx
index d55613a33b..51e07eb0e4 100644
--- a/src/libpsc/psc_push_particles/inc_curr_zigzag.cxx
+++ b/src/libpsc/psc_push_particles/inc_curr_zigzag.cxx
@@ -1,6 +1,10 @@
 
 #pragma once
 
+// T. Umeda, Y. Omura, T. Tominaga, and H. Matsumoto, "A new charge conservation
+// method in electromagnetic particle-in-cell simulations", Computer Physics
+// Communications 156 (2003) 73
+
 // ======================================================================
 
 template <typename Order, typename Dim, typename _fields_t>
diff --git a/src/libpsc/psc_push_particles/push_particles_1vb.hxx b/src/libpsc/psc_push_particles/push_particles_1vb.hxx
index 18c8ff1803..d2c0803003 100644
--- a/src/libpsc/psc_push_particles/push_particles_1vb.hxx
+++ b/src/libpsc/psc_push_particles/push_particles_1vb.hxx
@@ -1,4 +1,7 @@
 
+// J. Villasenor and O. Buneman, "Rigorous charge conservation for local
+// electromagnetic field solvers", Computer Physics Communications 69 (1992) 306
+
 // ======================================================================
 // PushParticlesVb
 
diff --git a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx
index 0b985f5cf8..58285a867e 100644
--- a/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx
+++ b/src/libpsc/psc_push_particles/push_particles_esirkepov.hxx
@@ -1,6 +1,10 @@
 
 #pragma once
 
+// T. Zh. Esirkepov, "Exact charge conservation scheme for particle-in-cell
+// simulation with an arbitrary form-factor", Computer Physics Communications
+// 135 (2001) 144
+
 #include "pushp_current_esirkepov.hxx"
 #include "../libpsc/psc_checks/checks_impl.hxx"
 
diff --git a/src/libpsc/psc_push_particles/pushp_current_esirkepov.hxx b/src/libpsc/psc_push_particles/pushp_current_esirkepov.hxx
index 0a25bac894..8e611ebc07 100644
--- a/src/libpsc/psc_push_particles/pushp_current_esirkepov.hxx
+++ b/src/libpsc/psc_push_particles/pushp_current_esirkepov.hxx
@@ -1,6 +1,10 @@
 
 #pragma once
 
+// T. Zh. Esirkepov, "Exact charge conservation scheme for particle-in-cell
+// simulation with an arbitrary form-factor", Computer Physics Communications
+// 135 (2001) 144
+
 #include "inc_defs.h"
 #include "interpolate.hxx"
 
diff --git a/src/libpsc/tests/test_current_deposition.cxx b/src/libpsc/tests/test_current_deposition.cxx
index ccd0f89816..b17f11d591 100644
--- a/src/libpsc/tests/test_current_deposition.cxx
+++ b/src/libpsc/tests/test_current_deposition.cxx
@@ -51,7 +51,7 @@ template <typename T>
 struct CurrentDepositionTest : ::testing::Test
 {
   using real_t = typename T::real_t;
-  using real3_t = Vec3<real_t>;
+  using Real3 = Vec3<real_t>;
   using dim_t = typename T::dim_t;
   using Current = typename T::Current;
 
@@ -74,12 +74,12 @@ struct CurrentDepositionTest : ::testing::Test
     grid_.reset(
       new Grid_t{grid_domain, grid_bc, kinds, grid_norm, dt, 1, this->ibn});
 
-    ldims_ = gt::shape_type<3>(grid().ldims);
+    ldims_ = Int3(grid().ldims);
   }
 
   const Grid_t& grid() const { return *grid_; }
 
-  auto calc_current(real3_t xm, real3_t xp, real3_t vxi)
+  auto calc_current(Real3 xm, Real3 xp, Real3 vxi)
   {
     using fields_view_t = typename T::fields_view_t;
     using fields_t = typename T::fields_t;
@@ -101,11 +101,11 @@ struct CurrentDepositionTest : ::testing::Test
   }
 
   template <typename F>
-  void check_continuity(real3_t xm, real3_t xp, const F& flds)
+  void check_continuity(Real3 xm, Real3 xp, const F& flds)
   {
     const real_t eps = 2. * std::numeric_limits<real_t>::epsilon();
 
-    auto ib = gt::shape(0, 0, 0);
+    Int3 ib{0, 0, 0};
     auto rho_m = gt::zeros<real_t>(ldims_);
     auto rho_p = gt::zeros<real_t>(ldims_);
     psc::deposit::norm::nc<dim_t>(rho_m, ib, xm, 1.);
@@ -125,7 +125,7 @@ struct CurrentDepositionTest : ::testing::Test
   }
 
   template <typename F>
-  void test_current(real3_t xm, real3_t xp, real3_t vxi, const F& jxi_ref,
+  void test_current(Real3 xm, Real3 xp, Real3 vxi, const F& jxi_ref,
                     const F& jyi_ref, const F& jzi_ref)
   {
     if (!std::is_same<dim_t, dim_yz>::value) {
@@ -143,7 +143,7 @@ struct CurrentDepositionTest : ::testing::Test
       << flds.view(0, _all, _all, JZI);
   }
 
-  void test_current(real3_t xm, real3_t xp, real3_t vxi)
+  void test_current(Real3 xm, Real3 xp, Real3 vxi)
   {
     auto flds = calc_current(xm, xp, vxi);
     check_continuity(xm, xp, flds);
@@ -151,7 +151,7 @@ struct CurrentDepositionTest : ::testing::Test
 
   Int3 ibn = {0, 0, 0};
   std::unique_ptr<Grid_t> grid_;
-  gt::shape_type<3> ldims_;
+  Int3 ldims_;
 };
 
 template <typename R, typename D,
@@ -180,10 +180,10 @@ TYPED_TEST(CurrentDepositionTest, CurrentNotMoving)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {.5, 1., 1.}, xp = {.5, 1., 1.};
-  real3_t vxi = {0., 0., 0.};
+  Real3 xm = {.5, 1., 1.}, xp = {.5, 1., 1.};
+  Real3 vxi = {0., 0., 0.};
   // clang-format off
   gt::gtensor<real_t, 2> jxi_ref = {{0., 0., 0., 0.},
                                     {0., 0., 0., 0.},
@@ -205,10 +205,10 @@ TYPED_TEST(CurrentDepositionTest, CurrentY)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {.5, 1., 1.}, xp = {.5, 1.2, 1.};
-  real3_t vxi = {0., .2, 0.};
+  Real3 xm = {.5, 1., 1.}, xp = {.5, 1.2, 1.};
+  Real3 vxi = {0., .2, 0.};
   // clang-format off
   gt::gtensor<real_t, 2> jxi_ref = {{0., 0., 0., 0.},
                                     {0., 0., 0., 0.},
@@ -230,10 +230,10 @@ TYPED_TEST(CurrentDepositionTest, CurrentYShift)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {.5, 1., 1.7}, xp = {.5, 1.2, 1.7};
-  real3_t vxi = {0., .2, 0.};
+  Real3 xm = {.5, 1., 1.7}, xp = {.5, 1.2, 1.7};
+  Real3 vxi = {0., .2, 0.};
   // clang-format off
   gt::gtensor<real_t, 2> jxi_ref = {{0., 0., 0., 0.},
                                     {0., 0., 0., 0.},
@@ -255,10 +255,10 @@ TYPED_TEST(CurrentDepositionTest, CurrentZ)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {.5, 1., 1.3}, xp = {.5, 1., 1.6};
-  real3_t vxi = {0., 0., 3.};
+  Real3 xm = {.5, 1., 1.3}, xp = {.5, 1., 1.6};
+  Real3 vxi = {0., 0., 3.};
   // clang-format off
   gt::gtensor<real_t, 2> jxi_ref = {{0., 0., 0., 0.},
                                     {0., 0., 0., 0.},
@@ -280,11 +280,11 @@ TYPED_TEST(CurrentDepositionTest, CurrentYCross)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {.5, 1.9, 1.}, xp = {.5, 2.1, 1.};
+  Real3 xm = {.5, 1.9, 1.}, xp = {.5, 2.1, 1.};
   Int3 lf = {0, 1, 1}, lg = {0, 1, 1};
-  real3_t vxi = {0., .2, 0.};
+  Real3 vxi = {0., .2, 0.};
   // clang-format off
   gt::gtensor<real_t, 2> jxi_ref = {{0., 0., 0., 0.},
                                     {0., 0., 0., 0.},
@@ -306,11 +306,11 @@ TYPED_TEST(CurrentDepositionTest, CurrentYCrossShift)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {.5, 1.9, 1.3}, xp = {.5, 2.1, 1.3};
+  Real3 xm = {.5, 1.9, 1.3}, xp = {.5, 2.1, 1.3};
   Int3 lf = {0, 1, 1}, lg = {0, 1, 1};
-  real3_t vxi = {0., .2, 0.};
+  Real3 vxi = {0., .2, 0.};
   // clang-format off
   gt::gtensor<real_t, 2> jxi_ref = {{0., 0., 0., 0.},
                                     {0., 0., 0., 0.},
@@ -332,11 +332,11 @@ TYPED_TEST(CurrentDepositionTest, CurrentYZ)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {.5, 1.2, 1.1}, xp = {.5, 1.4, 1.4};
+  Real3 xm = {.5, 1.2, 1.1}, xp = {.5, 1.4, 1.4};
   Int3 lf = {0, 1, 1}, lg = {0, 1, 1};
-  real3_t vxi = {0., .2, .3};
+  Real3 vxi = {0., .2, .3};
   // clang-format off
   gt::gtensor<real_t, 2> jxi_ref = {{0., 0., 0., 0.},
                                     {0., 0., 0., 0.},
@@ -358,11 +358,11 @@ TYPED_TEST(CurrentDepositionTest, CurrentYZCrossShift)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {.5, 1.9, 1.3}, xp = {.5, 2.1, 1.4};
+  Real3 xm = {.5, 1.9, 1.3}, xp = {.5, 2.1, 1.4};
   Int3 lf = {0, 1, 1}, lg = {0, 1, 1};
-  real3_t vxi = {0., .2, 0.};
+  Real3 vxi = {0., .2, 0.};
   // clang-format off
   gt::gtensor<real_t, 2> jxi_ref = {{0., 0., 0., 0.},
                                     {0., 0., 0., 0.},
@@ -384,10 +384,10 @@ TYPED_TEST(CurrentDepositionTest, CurrentYZCrossYZ)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {.5, 1.9, 1.6}, xp = {.5, 2.1, 2.2};
-  real3_t vxi = {0., .2, 0.};
+  Real3 xm = {.5, 1.9, 1.6}, xp = {.5, 2.1, 2.2};
+  Real3 vxi = {0., .2, 0.};
   if (std::is_same<typename TypeParam::Current,
                    Current1vbSplit<opt_order_1st, dim_yz,
                                    typename TypeParam::fields_t>>::value) {
@@ -415,10 +415,10 @@ TYPED_TEST(CurrentDepositionTest, CurrentX)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {1.3, 1., 1.}, xp = {1.6, 1., 1.};
-  real3_t vxi = {.3, 0., 0.};
+  Real3 xm = {1.3, 1., 1.}, xp = {1.6, 1., 1.};
+  Real3 vxi = {.3, 0., 0.};
   this->test_current(xm, xp, vxi);
 }
 
@@ -426,10 +426,10 @@ TYPED_TEST(CurrentDepositionTest, CurrentXY)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {1.1, 1.2, 1.3}, xp = {1.4, 1.8, 1.3};
-  real3_t vxi = {.3, .6, 0.};
+  Real3 xm = {1.1, 1.2, 1.3}, xp = {1.4, 1.8, 1.3};
+  Real3 vxi = {.3, .6, 0.};
   this->test_current(xm, xp, vxi);
 }
 
@@ -437,10 +437,10 @@ TYPED_TEST(CurrentDepositionTest, CurrentXYZ)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {1.1, 1.2, 1.3}, xp = {1.4, 1.6, 1.8};
-  real3_t vxi = {.3, .4, .5};
+  Real3 xm = {1.1, 1.2, 1.3}, xp = {1.4, 1.6, 1.8};
+  Real3 vxi = {.3, .4, .5};
   this->test_current(xm, xp, vxi);
 }
 
@@ -448,10 +448,10 @@ TYPED_TEST(CurrentDepositionTest, CurrentXYZCrossXYZ)
 {
   using self_type = CurrentDepositionTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
 
-  real3_t xm = {1.9, 1.8, 1.7}, xp = {2.2, 2.4, 2.6};
-  real3_t vxi = {.3, .6, .9};
+  Real3 xm = {1.9, 1.8, 1.7}, xp = {2.2, 2.4, 2.6};
+  Real3 vxi = {.3, .6, .9};
   this->test_current(xm, xp, vxi);
 }
 
diff --git a/src/libpsc/tests/test_deposit.cxx b/src/libpsc/tests/test_deposit.cxx
index e682ef5664..66debee107 100644
--- a/src/libpsc/tests/test_deposit.cxx
+++ b/src/libpsc/tests/test_deposit.cxx
@@ -16,7 +16,7 @@ template <typename T>
 struct DepositTest : ::testing::Test
 {
   using real_t = typename T::real_t;
-  using real3_t = gt::sarray<real_t, 3>;
+  using Real3 = Vec3<real_t>;
   using dim_t = typename T::dim_t;
 
   DepositTest() : ldims_{4, 4, 4}
@@ -26,7 +26,7 @@ struct DepositTest : ::testing::Test
     }
   }
 
-  gt::shape_type<3> ldims_;
+  Int3 ldims_;
   const real_t eps = std::numeric_limits<real_t>::epsilon();
 };
 
@@ -47,12 +47,12 @@ TYPED_TEST(DepositTest, ChargeCenter)
 {
   using self_type = DepositTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
   using dim_t = typename self_type::dim_t;
 
-  real3_t x = {1.5, 1.5, 1.5};
+  Real3 x = {1.5, 1.5, 1.5};
   real_t val = .1;
-  auto ibn = gt::shape(0, 0, 0);
+  Int3 ibn = {0, 0, 0};
   auto rho_ref = gt::zeros<real_t>(this->ldims_ + 2 * ibn);
   if (std::is_same<dim_t, dim_xyz>::value) {
     rho_ref.view(_s(1, 3), _s(1, 3), _s(1, 3)) = val / 8.f;
@@ -70,12 +70,12 @@ TYPED_TEST(DepositTest, ChargeLowerLeft)
 {
   using self_type = DepositTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
   using dim_t = typename self_type::dim_t;
 
-  real3_t x = {1., 1., 1.};
+  Real3 x = {1., 1., 1.};
   real_t val = .1;
-  auto ibn = gt::shape(0, 0, 0);
+  Int3 ibn = {0, 0, 0};
   auto rho_ref = gt::zeros<real_t>(this->ldims_ + 2 * ibn);
   if (std::is_same<dim_t, dim_xyz>::value) {
     rho_ref(1, 1, 1) = val;
@@ -93,12 +93,12 @@ TYPED_TEST(DepositTest, ChargeCenterWithBnd)
 {
   using self_type = DepositTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
   using dim_t = typename self_type::dim_t;
 
-  real3_t x = {.5, .5, .5};
+  Real3 x = {.5, .5, .5};
   real_t val = 1.;
-  auto ibn = gt::shape(1, 1, 1);
+  Int3 ibn = {1, 1, 1};
   auto rho_ref = gt::zeros<real_t>(this->ldims_ + 2 * ibn);
   if (std::is_same<dim_t, dim_xyz>::value) {
     rho_ref.view(_s(1, 3), _s(1, 3), _s(1, 3)) = val / 8.f;
@@ -116,12 +116,12 @@ TYPED_TEST(DepositTest, ChargeCenterWithBndCc)
 {
   using self_type = DepositTest<TypeParam>;
   using real_t = typename self_type::real_t;
-  using real3_t = typename self_type::real3_t;
+  using Real3 = typename self_type::Real3;
   using dim_t = typename self_type::dim_t;
 
-  real3_t x = {.5, .5, .5};
+  Real3 x = {.5, .5, .5};
   real_t val = 1.;
-  auto ibn = gt::shape(1, 1, 1);
+  Int3 ibn = {1, 1, 1};
   auto rho_ref = gt::zeros<real_t>(this->ldims_ + 2 * ibn);
   if (std::is_same<dim_t, dim_xyz>::value) {
     rho_ref(1, 1, 1) = val;