diff --git a/src/Expression/IpplExpressions.h b/src/Expression/IpplExpressions.h
index 632238eb2..fff85df4a 100644
--- a/src/Expression/IpplExpressions.h
+++ b/src/Expression/IpplExpressions.h
@@ -5,6 +5,8 @@
 #ifndef IPPL_EXPRESSIONS_H
 #define IPPL_EXPRESSIONS_H
 
+#include <Kokkos_Core.hpp>
+
 #include <type_traits>
 
 namespace ippl {
@@ -41,7 +43,7 @@ namespace ippl {
          * of the BareField class.
          */
         template <typename E, size_t N = sizeof(E)>
-        struct CapturedExpression {
+        struct alignas(E) CapturedExpression {
             constexpr static unsigned dim = E::dim;
 
             template <typename... Args>
@@ -50,7 +52,7 @@ namespace ippl {
                 return reinterpret_cast<const E&>(*this)(args...);
             }
 
-            char buffer[N];
+            alignas(E) char buffer[N];
         };
 
         /*!
@@ -94,6 +96,22 @@ namespace ippl {
         template <typename T>
         struct isExpression<Scalar<T>> : std::true_type {};
 
+        template <typename E, typename = void>
+        struct ExecutionSpaceOf {
+            using type = Kokkos::DefaultExecutionSpace;
+        };
+
+        template <typename E>
+        struct ExecutionSpaceOf<E, std::void_t<typename E::execution_space>> {
+            using type = typename E::execution_space;
+        };
+
+        template <typename E1, typename E2>
+        struct BinaryExecutionSpace {
+            using type = std::conditional_t<(E1::dim != 0), typename ExecutionSpaceOf<E1>::type,
+                                            typename ExecutionSpaceOf<E2>::type>;
+        };
+
     }  // namespace detail
 }  // namespace ippl
 
diff --git a/src/Expression/IpplOperations.h b/src/Expression/IpplOperations.h
index 1e02a6b8a..a63d07887 100644
--- a/src/Expression/IpplOperations.h
+++ b/src/Expression/IpplOperations.h
@@ -71,6 +71,7 @@ namespace ippl {
     struct fun : public detail::Expression<fun<E>, sizeof(E)> {                \
         constexpr static unsigned dim = E::dim;                                \
         using value_type              = typename E::value_type;                \
+        using execution_space         = typename detail::ExecutionSpaceOf<E>::type; \
                                                                                \
         KOKKOS_FUNCTION                                                        \
         fun(const E& u)                                                        \
@@ -132,6 +133,7 @@ namespace ippl {
     struct fun : public detail::Expression<fun<E1, E2>, sizeof(E1) + sizeof(E2)> {             \
         constexpr static unsigned dim = std::max(E1::dim, E2::dim);                            \
         using value_type              = typename E1::value_type;                               \
+        using execution_space         = typename detail::BinaryExecutionSpace<E1, E2>::type;   \
                                                                                                \
         KOKKOS_FUNCTION                                                                        \
         fun(const E1& u, const E2& v)                                                          \
diff --git a/src/Field/BareField.h b/src/Field/BareField.h
index 33908ef53..ff2c014d0 100644
--- a/src/Field/BareField.h
+++ b/src/Field/BareField.h
@@ -20,6 +20,7 @@
 
 #include "Field/HaloCells.h"
 #include "FieldLayout/FieldLayout.h"
+#include "Index/SIndex.h"
 
 namespace ippl {
     class Index;
@@ -166,6 +167,14 @@ namespace ippl {
             return dview_m(args...);
         }
 
+        auto operator[](const Domain_t& domain);
+
+        auto operator[](const Index& index);
+
+        auto operator[](int index);
+
+        auto operator[](const SIndex<Dim>& sindex);
+
         view_type& getView() { return dview_m; }
 
         const view_type& getView() const { return dview_m; }
diff --git a/src/Field/BareField.hpp b/src/Field/BareField.hpp
index efcddb5f1..44a3f1ff9 100644
--- a/src/Field/BareField.hpp
+++ b/src/Field/BareField.hpp
@@ -14,6 +14,8 @@
 
 #include "Utility/Inform.h"
 #include "Utility/IpplInfo.h"
+#include "Field/IndexedBareField.h"
+#include "Field/SparseIndexedBareField.h"
 namespace Kokkos {
     template <typename T, unsigned Dim>
     struct reduction_identity<ippl::Vector<T, Dim>> {
@@ -139,6 +141,29 @@ namespace ippl {
         Kokkos::resize(dview_m, args...);
     }
 
+    template <typename T, unsigned Dim, class... ViewArgs>
+    auto BareField<T, Dim, ViewArgs...>::operator[](const Domain_t& domain) {
+        return IndexedBareField<BareField<T, Dim, ViewArgs...>>(*this, domain);
+    }
+
+    template <typename T, unsigned Dim, class... ViewArgs>
+    auto BareField<T, Dim, ViewArgs...>::operator[](const Index& index) {
+        Domain_t domain = getDomain();
+        domain[0] = index;
+        return IndexedBareField<BareField<T, Dim, ViewArgs...>, 1>(*this, domain);
+    }
+
+    template <typename T, unsigned Dim, class... ViewArgs>
+    auto BareField<T, Dim, ViewArgs...>::operator[](int index) {
+        static_assert(Dim == 1, "Integer overload is only available for one-dimensional fields.");
+        return operator[](Index(index, index));
+    }
+
+    template <typename T, unsigned Dim, class... ViewArgs>
+    auto BareField<T, Dim, ViewArgs...>::operator[](const SIndex<Dim>& sindex) {
+        return SparseIndexedBareField<BareField<T, Dim, ViewArgs...>>(*this, sindex);
+    }
+
     template <typename T, unsigned Dim, class... ViewArgs>
     void BareField<T, Dim, ViewArgs...>::fillHalo() {
         if (layout_m->comm.size() > 1) {
diff --git a/src/Field/IndexedBareField.h b/src/Field/IndexedBareField.h
new file mode 100644
index 000000000..67e74abe5
--- /dev/null
+++ b/src/Field/IndexedBareField.h
@@ -0,0 +1,173 @@
+//
+// Class IndexedBareField
+//   Lightweight indexed view into a BareField.
+//
+#ifndef IPPL_INDEXED_BARE_FIELD_H
+#define IPPL_INDEXED_BARE_FIELD_H
+
+#include <Kokkos_Core.hpp>
+
+#include "Expression/IpplExpressions.h"
+#include "Index/NDIndex.h"
+#include "Utility/ParallelDispatch.h"
+
+namespace ippl {
+
+    template <typename Field, unsigned Brackets = Field::dim>
+    class IndexedBareField
+        : public detail::Expression<IndexedBareField<Field, Brackets>,
+                                    sizeof(Field) + sizeof(NDIndex<Field::dim>)> {
+    public:
+        using field_type      = Field;
+        using value_type      = typename field_type::value_type;
+        using view_type       = typename field_type::view_type;
+        using execution_space = typename field_type::execution_space;
+        using Domain_t        = NDIndex<field_type::dim>;
+
+        constexpr static unsigned dim = field_type::dim;
+        constexpr static size_t expression_size = sizeof(Field) + sizeof(Domain_t);
+
+        IndexedBareField(field_type& field, const Domain_t& domain)
+            : view_m(field.getView())
+            , owned_m(field.getOwned())
+            , domain_m(domain)
+            , nghost_m(field.getNghost()) {}
+
+        IndexedBareField(const IndexedBareField&) = default;
+
+        IndexedBareField(const view_type& view, const Domain_t& owned, const Domain_t& domain,
+                         int nghost)
+            : view_m(view)
+            , owned_m(owned)
+            , domain_m(domain)
+            , nghost_m(nghost) {}
+
+        auto operator[](const Index& index) const {
+            static_assert(Brackets < dim, "Too many Index brackets for field dimension.");
+            Domain_t domain = domain_m;
+            domain[Brackets] = index;
+            return IndexedBareField<Field, Brackets + 1>(view_m, owned_m, domain, nghost_m);
+        }
+
+        auto operator[](int index) const {
+            return operator[](Index(index, index));
+        }
+
+        template <typename... Args>
+        KOKKOS_INLINE_FUNCTION value_type operator()(Args... args) const {
+            static_assert(Brackets == dim, "IndexedBareField expression requires all dimensions.");
+            static_assert(sizeof...(Args) == dim);
+            typename RangePolicy<dim, execution_space>::index_array_type rel{args...};
+            return apply(view_m, relativeToView(rel));
+        }
+
+        IndexedBareField& operator=(value_type value) {
+            static_assert(Brackets == dim, "IndexedBareField assignment requires all dimensions.");
+            assign(ValueAssign{value});
+            return *this;
+        }
+
+        template <typename E, size_t N>
+        IndexedBareField& operator=(const detail::Expression<E, N>& expr) {
+            static_assert(Brackets == dim, "IndexedBareField assignment requires all dimensions.");
+            using capture_type = detail::CapturedExpression<E, N>;
+            capture_type expr_ = reinterpret_cast<const capture_type&>(expr);
+            assign(ExpressionAssign<E, N>{expr_});
+            return *this;
+        }
+
+        IndexedBareField& operator=(const IndexedBareField& rhs) {
+            static_assert(Brackets == dim, "IndexedBareField assignment requires all dimensions.");
+            const detail::Expression<IndexedBareField<Field, Brackets>, expression_size>& expr = rhs;
+            return operator=(expr);
+        }
+
+        const Domain_t& getDomain() const { return domain_m; }
+
+    private:
+        view_type view_m;
+        Domain_t owned_m;
+        Domain_t domain_m;
+        int nghost_m;
+
+    public:
+        struct ValueAssign {
+            value_type value;
+
+            template <typename Coords>
+            KOKKOS_INLINE_FUNCTION value_type operator()(const Coords&) const {
+                return value;
+            }
+        };
+
+        template <typename E, size_t N>
+        struct ExpressionAssign {
+            detail::CapturedExpression<E, N> expr;
+
+            template <typename Coords>
+            KOKKOS_INLINE_FUNCTION auto operator()(const Coords& rel) const {
+                return apply(expr, rel);
+            }
+        };
+
+    private:
+        template <typename Coords>
+        KOKKOS_INLINE_FUNCTION auto relativeToView(const Coords& rel) const {
+            typename RangePolicy<dim, execution_space>::index_array_type viewCoords;
+            for (unsigned d = 0; d < dim; ++d) {
+                const auto global = domain_m[d].first() + rel[d] * domain_m[d].stride();
+                viewCoords[d] = (global - owned_m[d].first()) / owned_m[d].stride() + nghost_m;
+            }
+            return viewCoords;
+        }
+
+        template <typename Coords>
+        KOKKOS_INLINE_FUNCTION auto globalToRelative(const Coords& global) const {
+            typename RangePolicy<dim, execution_space>::index_array_type rel;
+            for (unsigned d = 0; d < dim; ++d) {
+                rel[d] = (global[d] - domain_m[d].first()) / domain_m[d].stride();
+            }
+            return rel;
+        }
+
+        template <typename Coords>
+        KOKKOS_INLINE_FUNCTION auto globalToView(const Coords& global) const {
+            typename RangePolicy<dim, execution_space>::index_array_type viewCoords;
+            for (unsigned d = 0; d < dim; ++d) {
+                viewCoords[d] = (global[d] - owned_m[d].first()) / owned_m[d].stride() + nghost_m;
+            }
+            return viewCoords;
+        }
+
+    public:
+        template <typename Functor>
+        void assign(const Functor& functor) {
+            Domain_t local = domain_m.intersect(owned_m);
+            if (local.empty()) {
+                return;
+            }
+
+            using range_policy_type = RangePolicy<dim, execution_space>;
+            using index_type        = typename range_policy_type::index_type;
+            Kokkos::Array<index_type, dim> begin, end;
+            for (unsigned d = 0; d < dim; ++d) {
+                begin[d] = 0;
+                end[d]   = local[d].length();
+            }
+
+            ippl::parallel_for(
+                "IndexedBareField::operator=", createRangePolicy<dim, execution_space>(begin, end),
+                KOKKOS_CLASS_LAMBDA(const typename range_policy_type::index_array_type& args) {
+                    typename range_policy_type::index_array_type global;
+                    for (unsigned d = 0; d < dim; ++d) {
+                        global[d] = local[d].first() + args[d] * local[d].stride();
+                    }
+                    const auto rel = globalToRelative(global);
+                    apply(view_m, globalToView(global)) = functor(rel);
+                });
+        }
+    };
+
+}  // namespace ippl
+
+#endif
diff --git a/src/Field/SparseIndexedBareField.h b/src/Field/SparseIndexedBareField.h
new file mode 100644
index 000000000..5f102e3c7
--- /dev/null
+++ b/src/Field/SparseIndexedBareField.h
@@ -0,0 +1,273 @@
+//
+// Class SparseIndexedBareField
+//   Sparse point view into a BareField.
+//
+#ifndef IPPL_SPARSE_INDEXED_BARE_FIELD_H
+#define IPPL_SPARSE_INDEXED_BARE_FIELD_H
+
+#include <Kokkos_Core.hpp>
+
+#include "Expression/IpplExpressions.h"
+#include "Index/SIndex.h"
+#include "Utility/ParallelDispatch.h"
+
+namespace ippl {
+
+    template <typename Field>
+    class SparseIndexedBareField
+        : public detail::Expression<SparseIndexedBareField<Field>,
+                                    sizeof(Field) + sizeof(SIndex<Field::dim>)> {
+    public:
+        using field_type      = Field;
+        using value_type      = typename field_type::value_type;
+        using view_type       = typename field_type::view_type;
+        using execution_space = typename field_type::execution_space;
+        using sindex_type     = SIndex<field_type::dim>;
+        using point_type      = typename sindex_type::point_type;
+        using points_view_type = Kokkos::View<point_type*, execution_space>;
+
+        constexpr static unsigned dim = 1;
+        constexpr static unsigned field_dim = field_type::dim;
+        constexpr static size_t expression_size = sizeof(Field) + sizeof(sindex_type);
+
+        SparseIndexedBareField(field_type& field, const sindex_type& sindex)
+            : view_m(field.getView())
+            , owned_m(field.getOwned())
+            , points_m(sindex.template getDevicePoints<execution_space>())
+            , nghost_m(field.getNghost()) {}
+
+        SparseIndexedBareField(const SparseIndexedBareField&) = default;
+
+        template <typename Idx>
+        KOKKOS_INLINE_FUNCTION value_type operator()(Idx i) const {
+            return apply(view_m, pointToView(points_m(i)));
+        }
+
+        SparseIndexedBareField& operator=(value_type value) {
+            assign(ValueAssign{value});
+            return *this;
+        }
+
+        template <typename E, size_t N>
+        SparseIndexedBareField& operator=(const detail::Expression<E, N>& expr) {
+            using capture_type = detail::CapturedExpression<E, N>;
+            capture_type expr_ = reinterpret_cast<const capture_type&>(expr);
+            assign(ExpressionAssign<E, N>{expr_});
+            return *this;
+        }
+
+        SparseIndexedBareField& operator=(const SparseIndexedBareField& rhs) {
+            const detail::Expression<SparseIndexedBareField, expression_size>& expr = rhs;
+            return operator=(expr);
+        }
+
+        SparseIndexedBareField& operator+=(value_type value) {
+            update(PlusValueAssign{value});
+            return *this;
+        }
+
+        SparseIndexedBareField& operator-=(value_type value) {
+            update(MinusValueAssign{value});
+            return *this;
+        }
+
+        SparseIndexedBareField& operator*=(value_type value) {
+            update(MultiplyValueAssign{value});
+            return *this;
+        }
+
+        SparseIndexedBareField& operator/=(value_type value) {
+            update(DivideValueAssign{value});
+            return *this;
+        }
+
+        template <typename E, size_t N>
+        SparseIndexedBareField& operator+=(const detail::Expression<E, N>& expr) {
+            using capture_type = detail::CapturedExpression<E, N>;
+            capture_type expr_ = reinterpret_cast<const capture_type&>(expr);
+            update(PlusExpressionAssign<E, N>{expr_});
+            return *this;
+        }
+
+        template <typename E, size_t N>
+        SparseIndexedBareField& operator-=(const detail::Expression<E, N>& expr) {
+            using capture_type = detail::CapturedExpression<E, N>;
+            capture_type expr_ = reinterpret_cast<const capture_type&>(expr);
+            update(MinusExpressionAssign<E, N>{expr_});
+            return *this;
+        }
+
+        template <typename E, size_t N>
+        SparseIndexedBareField& operator*=(const detail::Expression<E, N>& expr) {
+            using capture_type = detail::CapturedExpression<E, N>;
+            capture_type expr_ = reinterpret_cast<const capture_type&>(expr);
+            update(MultiplyExpressionAssign<E, N>{expr_});
+            return *this;
+        }
+
+        template <typename E, size_t N>
+        SparseIndexedBareField& operator/=(const detail::Expression<E, N>& expr) {
+            using capture_type = detail::CapturedExpression<E, N>;
+            capture_type expr_ = reinterpret_cast<const capture_type&>(expr);
+            update(DivideExpressionAssign<E, N>{expr_});
+            return *this;
+        }
+
+    private:
+        view_type view_m;
+        NDIndex<field_dim> owned_m;
+        points_view_type points_m;
+        int nghost_m;
+
+    public:
+        struct ValueAssign {
+            value_type value;
+
+            template <typename Coords>
+            KOKKOS_INLINE_FUNCTION value_type operator()(const Coords&) const {
+                return value;
+            }
+        };
+
+        struct PlusValueAssign {
+            value_type value;
+
+            template <typename Current, typename Coords>
+            KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords&) const {
+                return current + value;
+            }
+        };
+
+        struct MinusValueAssign {
+            value_type value;
+
+            template <typename Current, typename Coords>
+            KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords&) const {
+                return current - value;
+            }
+        };
+
+        struct MultiplyValueAssign {
+            value_type value;
+
+            template <typename Current, typename Coords>
+            KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords&) const {
+                return current * value;
+            }
+        };
+
+        struct DivideValueAssign {
+            value_type value;
+
+            template <typename Current, typename Coords>
+            KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords&) const {
+                return current / value;
+            }
+        };
+
+        template <typename E, size_t N>
+        struct ExpressionAssign {
+            detail::CapturedExpression<E, N> expr;
+
+            template <typename Coords>
+            KOKKOS_INLINE_FUNCTION auto operator()(const Coords& rel) const {
+                return apply(expr, rel);
+            }
+        };
+
+        template <typename E, size_t N>
+        struct PlusExpressionAssign {
+            detail::CapturedExpression<E, N> expr;
+
+            template <typename Current, typename Coords>
+            KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords& rel) const {
+                return current + apply(expr, rel);
+            }
+        };
+
+        template <typename E, size_t N>
+        struct MinusExpressionAssign {
+            detail::CapturedExpression<E, N> expr;
+
+            template <typename Current, typename Coords>
+            KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords& rel) const {
+                return current - apply(expr, rel);
+            }
+        };
+
+        template <typename E, size_t N>
+        struct MultiplyExpressionAssign {
+            detail::CapturedExpression<E, N> expr;
+
+            template <typename Current, typename Coords>
+            KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords& rel) const {
+                return current * apply(expr, rel);
+            }
+        };
+
+        template <typename E, size_t N>
+        struct DivideExpressionAssign {
+            detail::CapturedExpression<E, N> expr;
+
+            template <typename Current, typename Coords>
+            KOKKOS_INLINE_FUNCTION auto operator()(const Current& current, const Coords& rel) const {
+                return current / apply(expr, rel);
+            }
+        };
+
+    private:
+        KOKKOS_INLINE_FUNCTION bool isLocal(const point_type& point) const {
+            for (unsigned d = 0; d < field_dim; ++d) {
+                if (!owned_m[d].contains(Index(point[d], point[d]))) {
+                    return false;
+                }
+            }
+            return true;
+        }
+
+        KOKKOS_INLINE_FUNCTION auto pointToView(const point_type& point) const {
+            typename RangePolicy<field_dim, execution_space>::index_array_type viewCoords;
+            for (unsigned d = 0; d < field_dim; ++d) {
+                viewCoords[d] = (point[d] - owned_m[d].first()) / owned_m[d].stride() + nghost_m;
+            }
+            return viewCoords;
+        }
+
+    public:
+        template <typename Functor>
+        void assign(const Functor& functor) {
+            ippl::parallel_for(
+                "SparseIndexedBareField::operator=", getRangePolicy(points_m),
+                KOKKOS_CLASS_LAMBDA(
+                    const typename RangePolicy<1, execution_space>::index_array_type& args) {
+                    const auto i = args[0];
+                    const point_type point = points_m(i);
+                    if (isLocal(point)) {
+                        typename RangePolicy<1, execution_space>::index_array_type rel;
+                        rel[0] = i;
+                        apply(view_m, pointToView(point)) = functor(rel);
+                    }
+                });
+        }
+
+        template <typename Functor>
+        void update(const Functor& functor) {
+            ippl::parallel_for(
+                "SparseIndexedBareField::operator op=", getRangePolicy(points_m),
+                KOKKOS_CLASS_LAMBDA(
+                    const typename RangePolicy<1, execution_space>::index_array_type& args) {
+                    const auto i = args[0];
+                    const point_type point = points_m(i);
+                    if (isLocal(point)) {
+                        typename RangePolicy<1, execution_space>::index_array_type rel;
+                        rel[0] = i;
+                        const auto viewCoords = pointToView(point);
+                        apply(view_m, viewCoords) = functor(apply(view_m, viewCoords), rel);
+                    }
+                });
+        }
+    };
+
+}  // namespace ippl
+
+#endif
diff --git a/src/Index/Index.h b/src/Index/Index.h
index b390f5f6e..13c725dd5 100644
--- a/src/Index/Index.h
+++ b/src/Index/Index.h
@@ -206,7 +206,7 @@ namespace ippl {
 
         // Test to see if there is any overlap between two Indexes.
         KOKKOS_INLINE_FUNCTION bool touches(const Index& a) const;
-        // Test to see if one contains another (endpoints only)
+        // Test to see if one contains all points of another.
         KOKKOS_INLINE_FUNCTION bool contains(const Index& a) const;
         // Split one into two.
         KOKKOS_INLINE_FUNCTION bool split(Index& l, Index& r) const;
diff --git a/src/Index/Index.hpp b/src/Index/Index.hpp
index 4c930c3fc..9c669ae88 100644
--- a/src/Index/Index.hpp
+++ b/src/Index/Index.hpp
@@ -152,11 +152,37 @@ namespace ippl {
     }
 
     KOKKOS_INLINE_FUNCTION bool Index::touches(const Index& a) const {
-        return (min() <= a.max()) && (max() >= a.min());
+        return !intersect(a).empty();
     }
 
     KOKKOS_INLINE_FUNCTION bool Index::contains(const Index& a) const {
-        return (min() <= a.min()) && (max() >= a.max());
+        if (a.empty()) {
+            return true;
+        }
+        const int thisMin = min();
+        const int thisMax = max();
+        const int thatMin = a.min();
+        const int thatMax = a.max();
+        int thisStride    = stride();
+        int thatStride    = a.stride();
+        if (thisStride < 0) {
+            thisStride = -thisStride;
+        }
+        if (thatStride < 0) {
+            thatStride = -thatStride;
+        }
+
+        const bool endpointsContained = (thisMin <= thatMin) && (thisMax >= thatMax);
+        if (!endpointsContained || thisStride == 1) {
+            return endpointsContained;
+        }
+
+        const bool endpointsOnStride =
+            ((thatMin - thisMin) % thisStride == 0) && ((thisMax - thatMax) % thisStride == 0);
+        if (a.length() == 1) {
+            return endpointsOnStride;
+        }
+        return endpointsOnStride && (thatStride % thisStride == 0);
     }
 
     KOKKOS_INLINE_FUNCTION bool Index::split(Index& l, Index& r) const {
@@ -261,7 +287,7 @@ namespace ippl {
     KOKKOS_INLINE_FUNCTION Index Index::grow(int ncells) const {
         Index index;
 
-        index.first_m  = this->first_m - ncells;
+        index.first_m  = this->first_m - ncells * this->stride_m;
         index.length_m = this->length_m + 2 * ncells;
         index.stride_m = this->stride_m;
 
diff --git a/src/Index/SIndex.h b/src/Index/SIndex.h
new file mode 100644
index 000000000..5abbefb05
--- /dev/null
+++ b/src/Index/SIndex.h
@@ -0,0 +1,468 @@
+//
+// Class SIndex
+//   Sparse set of field index points.
+//
+#ifndef IPPL_SINDEX_H
+#define IPPL_SINDEX_H
+
+#include <Kokkos_Core.hpp>
+
+#include <algorithm>
+#include <cstddef>
+#include <type_traits>
+#include <utility>
+#include <vector>
+
+#include "FieldLayout/FieldLayout.h"
+#include "Expression/IpplExpressions.h"
+#include "Index/NDIndex.h"
+#include "Index/SOffset.h"
+#include "Types/Vector.h"
+#include "Utility/ParallelDispatch.h"
+
+namespace ippl {
+
+    template <unsigned Dim>
+    class IndexedSIndex;
+
+    template <unsigned Dim>
+    class SIndex {
+    public:
+        using point_type = Vector<int, Dim>;
+
+        SIndex() = default;
+
+        explicit SIndex(FieldLayout<Dim>& layout)
+            : layout_m(&layout)
+            , domain_m(layout.getDomain()) {}
+
+        SIndex(const SIndex&) = default;
+
+        SIndex(const SIndex& sindex, const SOffset<Dim>& offset)
+            : layout_m(sindex.layout_m)
+            , domain_m(sindex.domain_m)
+            , points_m(sindex.points_m)
+            , offset_m(sindex.offset_m + offset) {}
+
+        void initialize(FieldLayout<Dim>& layout) {
+            layout_m = &layout;
+            domain_m = layout.getDomain();
+            points_m.clear();
+        }
+
+        bool needInitialize() const { return layout_m == nullptr; }
+
+        SIndex& operator=(const SIndex&) = default;
+
+        SIndex& operator=(const SOffset<Dim>& offset) {
+            clear();
+            offset_m = SOffset<Dim>();
+            addIndex(pointFromOffset(offset));
+            return *this;
+        }
+
+        SIndex& operator=(const NDIndex<Dim>& domain) {
+            clear();
+            addIndex(domain);
+            return *this;
+        }
+
+        template <typename E, size_t N>
+        SIndex& operator=(const detail::Expression<E, N>& expr) {
+            assignExpression(expr, domain_m, EvaluationMode::ViewCoordinates);
+            return *this;
+        }
+
+        SIndex& operator&=(const SIndex& rhs) {
+            eraseIf([&](const point_type& point) { return !rhs.hasIndex(point); });
+            domain_m = domain_m.intersect(rhs.domain_m);
+            return *this;
+        }
+
+        SIndex& operator&=(const SOffset<Dim>& offset) {
+            const point_type point = pointFromOffset(offset);
+            eraseIf([&](const point_type& existing) { return !samePoint(existing, point); });
+            return *this;
+        }
+
+        SIndex& operator&=(const NDIndex<Dim>& domain) {
+            eraseIf([&](const point_type& point) { return !contains(domain, effectivePoint(point)); });
+            domain_m = domain_m.intersect(domain);
+            return *this;
+        }
+
+        SIndex& operator|=(const SIndex& rhs) {
+            domain_m = boundingUnion(domain_m, rhs.domain_m);
+            for (const auto& point : rhs.points_m) {
+                addIndex(point);
+            }
+            return *this;
+        }
+
+        SIndex& operator|=(const SOffset<Dim>& offset) {
+            domain_m = boundingUnion(domain_m, singlePointDomain(pointFromOffset(offset)));
+            addIndex(pointFromOffset(offset));
+            return *this;
+        }
+
+        SIndex& operator|=(const NDIndex<Dim>& domain) {
+            domain_m = boundingUnion(domain_m, domain);
+            addIndex(domain);
+            return *this;
+        }
+
+        bool addIndex(const point_type& point) {
+            if (!contains(domain_m, point) || hasIndex(point)) {
+                return false;
+            }
+            points_m.push_back(point);
+            return true;
+        }
+
+        bool addIndex(const NDIndex<Dim>& point) {
+            if (point.empty()) {
+                return false;
+            }
+
+            bool addedAny = false;
+            point_type values;
+            addIndex(point, values, 0, addedAny);
+            return addedAny;
+        }
+
+        bool removeIndex(const point_type& point) {
+            const auto originalSize = points_m.size();
+            eraseIf([&](const point_type& existing) { return samePoint(existing, point); });
+            return points_m.size() != originalSize;
+        }
+
+        bool removeIndex(const SOffset<Dim>& offset) {
+            return removeIndex(pointFromOffset(offset));
+        }
+
+        bool removeIndex(const NDIndex<Dim>& point) {
+            for (unsigned d = 0; d < Dim; ++d) {
+                if (point[d].length() != 1) {
+                    return false;
+                }
+            }
+            return removeIndex(pointFromNDIndex(point));
+        }
+
+        void clear() { points_m.clear(); }
+
+        std::size_t size() const { return points_m.size(); }
+
+        bool empty() const { return points_m.empty(); }
+
+        const std::vector<point_type>& points() const { return points_m; }
+
+        const NDIndex<Dim>& getDomain() const { return domain_m; }
+
+        void setDomain(const NDIndex<Dim>& domain) { domain_m = domain; }
+
+        FieldLayout<Dim>& getFieldLayout() const { return *layout_m; }
+
+        SOffset<Dim>& getOffset() { return offset_m; }
+
+        const SOffset<Dim>& getOffset() const { return offset_m; }
+
+        bool hasIndex(const point_type& point) const {
+            for (const auto& existing : points_m) {
+                bool same = true;
+                for (unsigned d = 0; d < Dim; ++d) {
+                    same = same && existing[d] == point[d];
+                }
+                if (same) {
+                    return true;
+                }
+            }
+            return false;
+        }
+
+        bool hasIndex(const NDIndex<Dim>& point) const {
+            for (unsigned d = 0; d < Dim; ++d) {
+                if (point[d].length() != 1) {
+                    return false;
+                }
+            }
+            return hasIndex(pointFromNDIndex(point));
+        }
+
+        SIndex operator()(const SOffset<Dim>& offset) const { return SIndex(*this, offset); }
+
+        IndexedSIndex<Dim> operator[](const NDIndex<Dim>& domain) {
+            return IndexedSIndex<Dim>(*this, domain);
+        }
+
+        IndexedSIndex<Dim> operator[](const Index& index) {
+            NDIndex<Dim> domain = domain_m;
+            domain[0] = index;
+            return operator[](domain);
+        }
+
+        SIndex operator()(int i0) const {
+            static_assert(Dim == 1);
+            return operator()(SOffset<Dim>(i0));
+        }
+
+        SIndex operator()(int i0, int i1) const {
+            static_assert(Dim == 2);
+            return operator()(SOffset<Dim>(i0, i1));
+        }
+
+        SIndex operator()(int i0, int i1, int i2) const {
+            static_assert(Dim == 3);
+            return operator()(SOffset<Dim>(i0, i1, i2));
+        }
+
+        SIndex operator()(int i0, int i1, int i2, int i3) const {
+            static_assert(Dim == 4);
+            return operator()(SOffset<Dim>(i0, i1, i2, i3));
+        }
+
+        SIndex operator()(int i0, int i1, int i2, int i3, int i4) const {
+            static_assert(Dim == 5);
+            return operator()(SOffset<Dim>(i0, i1, i2, i3, i4));
+        }
+
+        SIndex operator()(int i0, int i1, int i2, int i3, int i4, int i5) const {
+            static_assert(Dim == 6);
+            return operator()(SOffset<Dim>(i0, i1, i2, i3, i4, i5));
+        }
+
+        template <typename ExecutionSpace>
+        Kokkos::View<point_type*, ExecutionSpace> getDevicePoints() const {
+            Kokkos::View<point_type*, ExecutionSpace> devicePoints("SIndex::points", points_m.size());
+            auto hostPoints = Kokkos::create_mirror_view(devicePoints);
+            for (std::size_t i = 0; i < points_m.size(); ++i) {
+                hostPoints(i) = effectivePoint(points_m[i]);
+            }
+            Kokkos::deep_copy(devicePoints, hostPoints);
+            return devicePoints;
+        }
+
+        friend SIndex operator+(const SIndex& sindex, const SOffset<Dim>& offset) {
+            return SIndex(sindex, offset);
+        }
+
+        friend SIndex operator+(const SOffset<Dim>& offset, const SIndex& sindex) {
+            return SIndex(sindex, offset);
+        }
+
+        friend SIndex operator-(const SIndex& sindex, const SOffset<Dim>& offset) {
+            return SIndex(sindex, -offset);
+        }
+
+    private:
+        FieldLayout<Dim>* layout_m = nullptr;
+        NDIndex<Dim> domain_m;
+        std::vector<point_type> points_m;
+        SOffset<Dim> offset_m;
+
+    public:
+        enum class EvaluationMode { ViewCoordinates, RelativeCoordinates };
+
+    private:
+        static bool contains(const NDIndex<Dim>& domain, const point_type& point) {
+            for (unsigned d = 0; d < Dim; ++d) {
+                if (!domain[d].contains(Index(point[d], point[d]))) {
+                    return false;
+                }
+            }
+            return true;
+        }
+
+        static bool samePoint(const point_type& lhs, const point_type& rhs) {
+            for (unsigned d = 0; d < Dim; ++d) {
+                if (lhs[d] != rhs[d]) {
+                    return false;
+                }
+            }
+            return true;
+        }
+
+        static point_type pointFromOffset(const SOffset<Dim>& offset) {
+            point_type point;
+            for (unsigned d = 0; d < Dim; ++d) {
+                point[d] = offset[d];
+            }
+            return point;
+        }
+
+        static point_type pointFromNDIndex(const NDIndex<Dim>& domain) {
+            point_type point;
+            for (unsigned d = 0; d < Dim; ++d) {
+                point[d] = domain[d].first();
+            }
+            return point;
+        }
+
+        point_type effectivePoint(const point_type& point) const {
+            point_type shifted = point;
+            for (unsigned d = 0; d < Dim; ++d) {
+                shifted[d] += offset_m[d];
+            }
+            return shifted;
+        }
+
+        template <typename Predicate>
+        void eraseIf(const Predicate& predicate) {
+            auto out = points_m.begin();
+            for (auto in = points_m.begin(); in != points_m.end(); ++in) {
+                if (!predicate(*in)) {
+                    *out = *in;
+                    ++out;
+                }
+            }
+            points_m.erase(out, points_m.end());
+        }
+
+        void addIndex(const NDIndex<Dim>& domain, point_type& values, unsigned d, bool& addedAny) {
+            if (d == Dim) {
+                addedAny = addIndex(values) || addedAny;
+                return;
+            }
+
+            for (int value = domain[d].first(); value <= domain[d].last();
+                 value += domain[d].stride()) {
+                values[d] = value;
+                addIndex(domain, values, d + 1, addedAny);
+            }
+        }
+
+        static NDIndex<Dim> boundingUnion(const NDIndex<Dim>& lhs, const NDIndex<Dim>& rhs) {
+            NDIndex<Dim> result;
+            for (unsigned d = 0; d < Dim; ++d) {
+                result[d] = Index(std::min(lhs[d].first(), rhs[d].first()),
+                                  std::max(lhs[d].last(), rhs[d].last()), lhs[d].stride());
+            }
+            return result;
+        }
+
+        static NDIndex<Dim> singlePointDomain(const point_type& point) {
+            NDIndex<Dim> result;
+            for (unsigned d = 0; d < Dim; ++d) {
+                result[d] = Index(point[d], point[d]);
+            }
+            return result;
+        }
+
+        template <typename Coords, typename IndexType>
+        KOKKOS_INLINE_FUNCTION static IndexType flatIndex(
+            const Coords& coords, const Kokkos::Array<IndexType, Dim>& extents) {
+            IndexType flat = 0;
+            for (unsigned d = 0; d < Dim; ++d) {
+                flat = flat * extents[d] + coords[d];
+            }
+            return flat;
+        }
+
+    public:
+        template <typename E, size_t N>
+        void assignExpression(const detail::Expression<E, N>& expr, const NDIndex<Dim>& domain,
+                              EvaluationMode mode) {
+            clear();
+            setDomain(domain);
+
+            NDIndex<Dim> local = domain.intersect(getFieldLayout().getLocalNDIndex());
+            if (local.empty()) {
+                return;
+            }
+
+            using execution_space   = typename detail::ExecutionSpaceOf<E>::type;
+            using range_policy_type = RangePolicy<Dim, execution_space>;
+            using index_type        = typename range_policy_type::index_type;
+            using capture_type      = detail::CapturedExpression<E, N>;
+
+            Kokkos::Array<index_type, Dim> begin, end, extents;
+            index_type size = 1;
+            for (unsigned d = 0; d < Dim; ++d) {
+                begin[d]   = 0;
+                end[d]     = local[d].length();
+                extents[d] = local[d].length();
+                size *= extents[d];
+            }
+
+            Kokkos::View<int*, execution_space> selected("SIndex::selected", size);
+            capture_type expr_ = reinterpret_cast<const capture_type&>(expr);
+            const int nghost   = 1;
+            auto localCopy     = local;
+            auto domainCopy    = domain;
+
+            ippl::parallel_for(
+                "SIndex::operator=", createRangePolicy<Dim, execution_space>(begin, end),
+                KOKKOS_LAMBDA(const typename range_policy_type::index_array_type& args) {
+                    typename range_policy_type::index_array_type exprCoords;
+                    for (unsigned d = 0; d < Dim; ++d) {
+                        const auto global = localCopy[d].first() + args[d] * localCopy[d].stride();
+                        if (mode == EvaluationMode::ViewCoordinates) {
+                            exprCoords[d] =
+                                (global - localCopy[d].first()) / localCopy[d].stride() + nghost;
+                        } else {
+                            exprCoords[d] =
+                                (global - domainCopy[d].first()) / domainCopy[d].stride();
+                        }
+                    }
+                    selected(flatIndex(args, extents)) = apply(expr_, exprCoords) ? 1 : 0;
+                });
+
+            auto hostSelected = Kokkos::create_mirror_view(selected);
+            Kokkos::deep_copy(hostSelected, selected);
+
+            std::size_t selectedCount = 0;
+            for (index_type flat = 0; flat < size; ++flat) {
+                selectedCount += hostSelected(flat) != 0;
+            }
+            points_m.reserve(selectedCount);
+
+            point_type global;
+            appendSelected(local, extents, hostSelected, global, 0, index_type{0});
+        }
+
+        template <typename HostView, typename IndexType>
+        void appendSelected(const NDIndex<Dim>& local, const Kokkos::Array<IndexType, Dim>& extents,
+                            const HostView& selected, point_type& global, unsigned d,
+                            IndexType flat) {
+            if (d == Dim) {
+                if (selected(flat)) {
+                    points_m.push_back(global);
+                }
+                return;
+            }
+
+            for (IndexType i = 0; i < extents[d]; ++i) {
+                global[d] = local[d].first() + static_cast<int>(i) * local[d].stride();
+                appendSelected(local, extents, selected, global, d + 1, flat * extents[d] + i);
+            }
+        }
+
+        friend class IndexedSIndex<Dim>;
+    };
+
+    template <unsigned Dim>
+    class IndexedSIndex {
+    public:
+        IndexedSIndex(SIndex<Dim>& sindex, const NDIndex<Dim>& domain)
+            : sindex_m(sindex)
+            , domain_m(domain) {}
+
+        template <typename E, size_t N>
+        IndexedSIndex& operator=(const detail::Expression<E, N>& expr) {
+            sindex_m.assignExpression(expr, domain_m, SIndex<Dim>::EvaluationMode::RelativeCoordinates);
+            return *this;
+        }
+
+    private:
+        SIndex<Dim>& sindex_m;
+        NDIndex<Dim> domain_m;
+    };
+
+    template <typename E, size_t N, typename T, typename = std::enable_if_t<std::is_scalar_v<T>>>
+    auto gt(const detail::Expression<E, N>& expr, const T& value) {
+        return expr > value;
+    }
+
+}  // namespace ippl
+
+#endif
diff --git a/src/Index/SOffset.h b/src/Index/SOffset.h
new file mode 100644
index 000000000..2a871ff65
--- /dev/null
+++ b/src/Index/SOffset.h
@@ -0,0 +1,133 @@
+//
+// Class SOffset
+//   Integer offset for sparse field index points.
+//
+#ifndef IPPL_SOFFSET_H
+#define IPPL_SOFFSET_H
+
+#include <ostream>
+#include <type_traits>
+
+#include "Index/NDIndex.h"
+#include "Types/Vector.h"
+
+namespace ippl {
+
+    template <unsigned Dim>
+    class SOffset {
+    public:
+        using point_type = Vector<int, Dim>;
+
+        SOffset() = default;
+
+        explicit SOffset(int value)
+            : values_m(value) {}
+
+        template <typename... Args,
+                  typename std::enable_if<sizeof...(Args) == Dim, bool>::type = true>
+        explicit SOffset(Args... args)
+            : values_m(static_cast<int>(args)...) {}
+
+        explicit SOffset(const point_type& values)
+            : values_m(values) {}
+
+        int& operator[](unsigned d) { return values_m[d]; }
+
+        int operator[](unsigned d) const { return values_m[d]; }
+
+        const point_type& values() const { return values_m; }
+
+        SOffset& operator+=(const SOffset& rhs) {
+            for (unsigned d = 0; d < Dim; ++d) {
+                values_m[d] += rhs[d];
+            }
+            return *this;
+        }
+
+        SOffset& operator-=(const SOffset& rhs) {
+            for (unsigned d = 0; d < Dim; ++d) {
+                values_m[d] -= rhs[d];
+            }
+            return *this;
+        }
+
+        bool operator==(const SOffset& rhs) const {
+            for (unsigned d = 0; d < Dim; ++d) {
+                if (values_m[d] != rhs[d]) {
+                    return false;
+                }
+            }
+            return true;
+        }
+
+        bool operator!=(const SOffset& rhs) const { return !(*this == rhs); }
+
+        bool inside(const NDIndex<Dim>& domain) const {
+            for (unsigned d = 0; d < Dim; ++d) {
+                if (!domain[d].contains(Index(values_m[d], values_m[d]))) {
+                    return false;
+                }
+            }
+            return true;
+        }
+
+    private:
+        point_type values_m;
+    };
+
+    template <unsigned Dim>
+    SOffset<Dim> operator+(SOffset<Dim> lhs, const SOffset<Dim>& rhs) {
+        lhs += rhs;
+        return lhs;
+    }
+
+    template <unsigned Dim>
+    SOffset<Dim> operator-(SOffset<Dim> lhs, const SOffset<Dim>& rhs) {
+        lhs -= rhs;
+        return lhs;
+    }
+
+    template <unsigned Dim>
+    SOffset<Dim> operator-(const SOffset<Dim>& value) {
+        SOffset<Dim> result;
+        for (unsigned d = 0; d < Dim; ++d) {
+            result[d] = -value[d];
+        }
+        return result;
+    }
+
+    template <unsigned Dim>
+    NDIndex<Dim> operator+(NDIndex<Dim> domain, const SOffset<Dim>& offset) {
+        for (unsigned d = 0; d < Dim; ++d) {
+            domain[d] = Index(domain[d].first() + offset[d], domain[d].last() + offset[d],
+                              domain[d].stride());
+        }
+        return domain;
+    }
+
+    template <unsigned Dim>
+    NDIndex<Dim> operator+(const SOffset<Dim>& offset, const NDIndex<Dim>& domain) {
+        return domain + offset;
+    }
+
+    template <unsigned Dim>
+    NDIndex<Dim> operator-(const NDIndex<Dim>& domain, const SOffset<Dim>& offset) {
+        return domain + (-offset);
+    }
+
+    template <unsigned Dim>
+    std::ostream& operator<<(std::ostream& out, const SOffset<Dim>& offset) {
+        out << "[";
+        for (unsigned d = 0; d < Dim; ++d) {
+            out << offset[d];
+            if (d + 1 < Dim) {
+                out << ",";
+            }
+        }
+        out << "]";
+        return out;
+    }
+
+}  // namespace ippl
+
+#endif
diff --git a/src/IpplCore.h b/src/IpplCore.h
index 4ee5d796c..51fffe3e2 100644
--- a/src/IpplCore.h
+++ b/src/IpplCore.h
@@ -18,6 +18,8 @@
 
 // #include "Index/Index.h"
 // #include "Index/NDIndex.h"
+#include "Index/SOffset.h"
+#include "Index/SIndex.h"
 
 #include "FieldLayout/FieldLayout.h"
 #include "FieldLayout/SubFieldLayout.h"
diff --git a/unit_tests/BareField/BareField.cpp b/unit_tests/BareField/BareField.cpp
index fe439dc6f..4bf8b0ecc 100644
--- a/unit_tests/BareField/BareField.cpp
+++ b/unit_tests/BareField/BareField.cpp
@@ -106,6 +106,530 @@ TYPED_TEST(BareFieldTest, DeepCopy) {
     });
 }
 
+TYPED_TEST(BareFieldTest, IndexedSubdomainScalarAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& field = this->field;
+    *field      = 0;
+
+    ippl::NDIndex<TestFixture::dim> subdomain = field->getDomain();
+    subdomain[0] =
+        ippl::Index(subdomain[0].first(), subdomain[0].last(), 2 * subdomain[0].stride());
+
+    (*field)[subdomain] = T(7);
+
+    auto mirror = field->getHostMirror();
+    Kokkos::deep_copy(mirror, field->getView());
+
+    const auto owned = field->getOwned();
+    const int nghost = field->getNghost();
+    nestedViewLoop(mirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        std::array<ippl::Index, TestFixture::dim> pointIndices;
+        const std::array<size_t, TestFixture::dim> viewCoords{static_cast<size_t>(args)...};
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            const int global = owned[d].first()
+                               + static_cast<int>(viewCoords[d] - nghost) * owned[d].stride();
+            pointIndices[d] = ippl::Index(global, global);
+        }
+        auto point = std::make_from_tuple<ippl::NDIndex<TestFixture::dim>>(pointIndices);
+        const T expected = subdomain.contains(point) ? T(7) : T(0);
+        assertEqual<T>(expected, mirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, IndexedSubdomainExpressionAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& lhs = this->field;
+    lhs->operator=(0);
+
+    typename TestFixture::field_type rhs(this->layout);
+    rhs = T(5);
+
+    ippl::NDIndex<TestFixture::dim> subdomain = lhs->getDomain();
+    subdomain[0] =
+        ippl::Index(subdomain[0].first(), subdomain[0].last(), 2 * subdomain[0].stride());
+
+    (*lhs)[subdomain] = rhs[subdomain] + T(2);
+
+    auto mirror = lhs->getHostMirror();
+    Kokkos::deep_copy(mirror, lhs->getView());
+
+    const auto owned = lhs->getOwned();
+    const int nghost = lhs->getNghost();
+    nestedViewLoop(mirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        std::array<ippl::Index, TestFixture::dim> pointIndices;
+        const std::array<size_t, TestFixture::dim> viewCoords{static_cast<size_t>(args)...};
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            const int global = owned[d].first()
+                               + static_cast<int>(viewCoords[d] - nghost) * owned[d].stride();
+            pointIndices[d] = ippl::Index(global, global);
+        }
+        auto point = std::make_from_tuple<ippl::NDIndex<TestFixture::dim>>(pointIndices);
+        const T expected = subdomain.contains(point) ? T(7) : T(0);
+        assertEqual<T>(expected, mirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, IndexedSubdomainDifferentStrideExpressionAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& lhs = this->field;
+    lhs->operator=(0);
+
+    typename TestFixture::field_type rhs(this->layout);
+    auto rhsView   = rhs.getView();
+    const auto lDom = rhs.getLayout().getLocalNDIndex();
+    Kokkos::parallel_for("Set rhs field", rhs.getFieldRangePolicy(), FieldVal<TypeParam>{rhsView, lDom});
+    Kokkos::fence();
+
+    ippl::NDIndex<TestFixture::dim> lhsSubdomain = lhs->getOwned();
+    lhsSubdomain[0] = ippl::Index(lhsSubdomain[0].first(), lhsSubdomain[0].last(),
+                                  2 * lhsSubdomain[0].stride());
+
+    ippl::NDIndex<TestFixture::dim> rhsSubdomain = rhs.getOwned();
+    const int rhsLast = rhsSubdomain[0].first()
+                        + (lhsSubdomain[0].length() - 1) * rhsSubdomain[0].stride();
+    rhsSubdomain[0] = ippl::Index(rhsSubdomain[0].first(), rhsLast, rhsSubdomain[0].stride());
+
+    ASSERT_EQ(lhsSubdomain[0].length(), rhsSubdomain[0].length());
+    ASSERT_NE(lhsSubdomain[0].stride(), rhsSubdomain[0].stride());
+
+    (*lhs)[lhsSubdomain] = rhs[rhsSubdomain] + T(2);
+
+    auto lhsMirror = lhs->getHostMirror();
+    auto rhsMirror = rhs.getHostMirror();
+    Kokkos::deep_copy(lhsMirror, lhs->getView());
+    Kokkos::deep_copy(rhsMirror, rhs.getView());
+
+    const auto owned = lhs->getOwned();
+    const int nghost = lhs->getNghost();
+    nestedViewLoop(lhsMirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        std::array<ippl::Index, TestFixture::dim> pointIndices;
+        typename ippl::RangePolicy<TestFixture::dim, typename TestFixture::exec_space>::index_array_type
+            rhsCoords;
+        const std::array<size_t, TestFixture::dim> viewCoords{static_cast<size_t>(args)...};
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            const int global = owned[d].first()
+                               + static_cast<int>(viewCoords[d] - nghost) * owned[d].stride();
+            pointIndices[d] = ippl::Index(global, global);
+
+            const int rel       = (global - lhsSubdomain[d].first()) / lhsSubdomain[d].stride();
+            const int rhsGlobal = rhsSubdomain[d].first() + rel * rhsSubdomain[d].stride();
+            rhsCoords[d] = (rhsGlobal - owned[d].first()) / owned[d].stride() + nghost;
+        }
+        auto point = std::make_from_tuple<ippl::NDIndex<TestFixture::dim>>(pointIndices);
+        const T expected = lhsSubdomain.contains(point) ? ippl::apply(rhsMirror, rhsCoords) + T(2)
+                                                        : T(0);
+        assertEqual<T>(expected, lhsMirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, ChainedIndexedSubdomainScalarAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& field = this->field;
+    *field      = 0;
+
+    ippl::NDIndex<TestFixture::dim> subdomain = field->getDomain();
+    subdomain[0] =
+        ippl::Index(subdomain[0].first(), subdomain[0].last(), 2 * subdomain[0].stride());
+
+    if constexpr (TestFixture::dim == 1) {
+        (*field)[subdomain[0]] = T(7);
+    } else if constexpr (TestFixture::dim == 2) {
+        (*field)[subdomain[0]][subdomain[1]] = T(7);
+    } else if constexpr (TestFixture::dim == 3) {
+        (*field)[subdomain[0]][subdomain[1]][subdomain[2]] = T(7);
+    } else if constexpr (TestFixture::dim == 4) {
+        (*field)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]] = T(7);
+    } else if constexpr (TestFixture::dim == 5) {
+        (*field)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]] = T(7);
+    } else if constexpr (TestFixture::dim == 6) {
+        (*field)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]]
+                [subdomain[5]] = T(7);
+    }
+
+    auto mirror = field->getHostMirror();
+    Kokkos::deep_copy(mirror, field->getView());
+
+    const auto owned = field->getOwned();
+    const int nghost = field->getNghost();
+    nestedViewLoop(mirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        std::array<ippl::Index, TestFixture::dim> pointIndices;
+        const std::array<size_t, TestFixture::dim> viewCoords{static_cast<size_t>(args)...};
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            const int global = owned[d].first()
+                               + static_cast<int>(viewCoords[d] - nghost) * owned[d].stride();
+            pointIndices[d] = ippl::Index(global, global);
+        }
+        auto point = std::make_from_tuple<ippl::NDIndex<TestFixture::dim>>(pointIndices);
+        const T expected = subdomain.contains(point) ? T(7) : T(0);
+        assertEqual<T>(expected, mirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, ChainedIndexedSubdomainExpressionAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& lhs = this->field;
+    lhs->operator=(0);
+
+    typename TestFixture::field_type rhs(this->layout);
+    rhs = T(5);
+
+    ippl::NDIndex<TestFixture::dim> subdomain = lhs->getDomain();
+    subdomain[0] =
+        ippl::Index(subdomain[0].first(), subdomain[0].last(), 2 * subdomain[0].stride());
+
+    if constexpr (TestFixture::dim == 1) {
+        (*lhs)[subdomain[0]] = rhs[subdomain[0]] + T(2);
+    } else if constexpr (TestFixture::dim == 2) {
+        (*lhs)[subdomain[0]][subdomain[1]] = rhs[subdomain[0]][subdomain[1]] + T(2);
+    } else if constexpr (TestFixture::dim == 3) {
+        (*lhs)[subdomain[0]][subdomain[1]][subdomain[2]] =
+            rhs[subdomain[0]][subdomain[1]][subdomain[2]] + T(2);
+    } else if constexpr (TestFixture::dim == 4) {
+        (*lhs)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]] =
+            rhs[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]] + T(2);
+    } else if constexpr (TestFixture::dim == 5) {
+        (*lhs)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]] =
+            rhs[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]] + T(2);
+    } else if constexpr (TestFixture::dim == 6) {
+        (*lhs)[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]]
+              [subdomain[5]] =
+                  rhs[subdomain[0]][subdomain[1]][subdomain[2]][subdomain[3]][subdomain[4]]
+                     [subdomain[5]]
+                  + T(2);
+    }
+
+    auto mirror = lhs->getHostMirror();
+    Kokkos::deep_copy(mirror, lhs->getView());
+
+    const auto owned = lhs->getOwned();
+    const int nghost = lhs->getNghost();
+    nestedViewLoop(mirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        std::array<ippl::Index, TestFixture::dim> pointIndices;
+        const std::array<size_t, TestFixture::dim> viewCoords{static_cast<size_t>(args)...};
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            const int global = owned[d].first()
+                               + static_cast<int>(viewCoords[d] - nghost) * owned[d].stride();
+            pointIndices[d] = ippl::Index(global, global);
+        }
+        auto point = std::make_from_tuple<ippl::NDIndex<TestFixture::dim>>(pointIndices);
+        const T expected = subdomain.contains(point) ? T(7) : T(0);
+        assertEqual<T>(expected, mirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, SparseIndexedSubdomainScalarAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& field = this->field;
+    *field      = 0;
+
+    ippl::SIndex<TestFixture::dim> sindex(this->layout);
+    const auto domain = field->getDomain();
+
+    for (int i = domain[0].first(); i <= domain[0].last(); i += 2 * domain[0].stride()) {
+        typename ippl::SIndex<TestFixture::dim>::point_type point;
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            point[d] = domain[d].first();
+        }
+        point[0] = i;
+        sindex.addIndex(point);
+    }
+
+    (*field)[sindex] = T(7);
+
+    auto mirror = field->getHostMirror();
+    Kokkos::deep_copy(mirror, field->getView());
+
+    const auto owned = field->getOwned();
+    const int nghost = field->getNghost();
+    nestedViewLoop(mirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        typename ippl::SIndex<TestFixture::dim>::point_type point;
+        const std::array<size_t, TestFixture::dim> viewCoords{static_cast<size_t>(args)...};
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            point[d] = owned[d].first()
+                       + static_cast<int>(viewCoords[d] - nghost) * owned[d].stride();
+        }
+        const T expected = sindex.hasIndex(point) ? T(7) : T(0);
+        assertEqual<T>(expected, mirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, SparseIndexedSubdomainExpressionAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& lhs = this->field;
+    lhs->operator=(0);
+
+    typename TestFixture::field_type rhs(this->layout);
+    rhs = T(5);
+
+    ippl::SIndex<TestFixture::dim> sindex(this->layout);
+    const auto domain = lhs->getDomain();
+
+    for (int i = domain[0].first(); i <= domain[0].last(); i += 2 * domain[0].stride()) {
+        typename ippl::SIndex<TestFixture::dim>::point_type point;
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            point[d] = domain[d].first();
+        }
+        point[0] = i;
+        sindex.addIndex(point);
+    }
+
+    (*lhs)[sindex] = rhs[sindex] + T(2);
+
+    auto mirror = lhs->getHostMirror();
+    Kokkos::deep_copy(mirror, lhs->getView());
+
+    const auto owned = lhs->getOwned();
+    const int nghost = lhs->getNghost();
+    nestedViewLoop(mirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        typename ippl::SIndex<TestFixture::dim>::point_type point;
+        const std::array<size_t, TestFixture::dim> viewCoords{static_cast<size_t>(args)...};
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            point[d] = owned[d].first()
+                       + static_cast<int>(viewCoords[d] - nghost) * owned[d].stride();
+        }
+        const T expected = sindex.hasIndex(point) ? T(7) : T(0);
+        assertEqual<T>(expected, mirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, SparseIndexedOffsetExpressionAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& lhs = this->field;
+    lhs->operator=(0);
+
+    typename TestFixture::field_type rhs(this->layout);
+    rhs = T(5);
+    rhs.fillHalo();
+
+    ippl::SIndex<TestFixture::dim> sindex(this->layout);
+    const auto domain = lhs->getDomain();
+    const int offsetValue = domain[0].stride();
+
+    for (int i = domain[0].first(); i <= domain[0].last() - offsetValue;
+         i += 2 * domain[0].stride()) {
+        typename ippl::SIndex<TestFixture::dim>::point_type point;
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            point[d] = domain[d].first();
+        }
+        point[0] = i;
+        sindex.addIndex(point);
+    }
+
+    ippl::SOffset<TestFixture::dim> offset;
+    offset[0] = offsetValue;
+
+    (*lhs)[sindex] = rhs[sindex + offset] + T(2);
+
+    auto mirror = lhs->getHostMirror();
+    Kokkos::deep_copy(mirror, lhs->getView());
+
+    const auto owned = lhs->getOwned();
+    const int nghost = lhs->getNghost();
+    nestedViewLoop(mirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        typename ippl::SIndex<TestFixture::dim>::point_type point;
+        const std::array<size_t, TestFixture::dim> viewCoords{static_cast<size_t>(args)...};
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            point[d] = owned[d].first()
+                       + static_cast<int>(viewCoords[d] - nghost) * owned[d].stride();
+        }
+        const T expected = sindex.hasIndex(point) ? T(7) : T(0);
+        assertEqual<T>(expected, mirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, SparseIndexedCallOffsetExpressionAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& lhs = this->field;
+    lhs->operator=(0);
+
+    typename TestFixture::field_type rhs(this->layout);
+    rhs = T(5);
+    rhs.fillHalo();
+
+    ippl::SIndex<TestFixture::dim> sindex(this->layout);
+    const auto domain = lhs->getDomain();
+    const int offsetValue = domain[0].stride();
+
+    for (int i = domain[0].first(); i <= domain[0].last() - offsetValue;
+         i += 2 * domain[0].stride()) {
+        typename ippl::SIndex<TestFixture::dim>::point_type point;
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            point[d] = domain[d].first();
+        }
+        point[0] = i;
+        sindex.addIndex(point);
+    }
+
+    ippl::SOffset<TestFixture::dim> offset;
+    offset[0] = offsetValue;
+
+    (*lhs)[sindex] = rhs[sindex(offset)] + T(2);
+
+    auto mirror = lhs->getHostMirror();
+    Kokkos::deep_copy(mirror, lhs->getView());
+
+    const auto owned = lhs->getOwned();
+    const int nghost = lhs->getNghost();
+    nestedViewLoop(mirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        typename ippl::SIndex<TestFixture::dim>::point_type point;
+        const std::array<size_t, TestFixture::dim> viewCoords{static_cast<size_t>(args)...};
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            point[d] = owned[d].first()
+                       + static_cast<int>(viewCoords[d] - nghost) * owned[d].stride();
+        }
+        const T expected = sindex.hasIndex(point) ? T(7) : T(0);
+        assertEqual<T>(expected, mirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, SparseIndexSetOperations) {
+    ippl::SIndex<TestFixture::dim> sindex(this->layout);
+    const auto domain = this->field->getDomain();
+
+    typename ippl::SIndex<TestFixture::dim>::point_type firstPoint;
+    typename ippl::SIndex<TestFixture::dim>::point_type secondPoint;
+    for (unsigned d = 0; d < TestFixture::dim; ++d) {
+        firstPoint[d]  = domain[d].first();
+        secondPoint[d] = domain[d].first();
+    }
+    secondPoint[0] = domain[0].first() + domain[0].stride();
+
+    EXPECT_TRUE(sindex.addIndex(firstPoint));
+    EXPECT_TRUE(sindex.addIndex(secondPoint));
+    EXPECT_FALSE(sindex.addIndex(firstPoint));
+
+    ippl::NDIndex<TestFixture::dim> subset = domain;
+    subset[0] = ippl::Index(firstPoint[0], firstPoint[0], domain[0].stride());
+    sindex &= subset;
+
+    EXPECT_TRUE(sindex.hasIndex(firstPoint));
+    EXPECT_FALSE(sindex.hasIndex(secondPoint));
+
+    ippl::SOffset<TestFixture::dim> secondOffset;
+    for (unsigned d = 0; d < TestFixture::dim; ++d) {
+        secondOffset[d] = secondPoint[d];
+    }
+    sindex |= secondOffset;
+
+    EXPECT_TRUE(sindex.hasIndex(secondPoint));
+    EXPECT_EQ(sindex.size(), 2u);
+}
+
+TYPED_TEST(BareFieldTest, SparseIndexPredicateAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& source = this->field;
+    auto view    = source->getView();
+    const auto lDom = source->getLayout().getLocalNDIndex();
+
+    Kokkos::parallel_for("Set field", source->getFieldRangePolicy(),
+                         FieldVal<TypeParam>{view, lDom});
+    Kokkos::fence();
+
+    typename TestFixture::field_type selected(this->layout);
+    selected = T(0);
+
+    ippl::SIndex<TestFixture::dim> sindex(this->layout);
+    sindex = *source > T(2);
+    selected[sindex] = (*source)[sindex];
+
+    auto sourceMirror = source->getHostMirror();
+    auto selectedMirror = selected.getHostMirror();
+    Kokkos::deep_copy(sourceMirror, source->getView());
+    Kokkos::deep_copy(selectedMirror, selected.getView());
+
+    const int nghost = source->getNghost();
+    nestedViewLoop(sourceMirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        const T value = sourceMirror(args...);
+        const T expected = value > T(2) ? value : T(0);
+        assertEqual<T>(expected, selectedMirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, SparseIndexSubsetPredicateAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& source = this->field;
+    auto view    = source->getView();
+    const auto lDom = source->getLayout().getLocalNDIndex();
+
+    Kokkos::parallel_for("Set field", source->getFieldRangePolicy(),
+                         FieldVal<TypeParam>{view, lDom});
+    Kokkos::fence();
+
+    typename TestFixture::field_type selected(this->layout);
+    selected = T(0);
+
+    ippl::NDIndex<TestFixture::dim> subset = source->getDomain();
+    subset[0] = ippl::Index(subset[0].first(), subset[0].last(), 2 * subset[0].stride());
+
+    ippl::SIndex<TestFixture::dim> sindex(this->layout);
+    sindex[subset] = (*source)[subset] > T(2);
+    selected[sindex] = (*source)[sindex];
+
+    auto sourceMirror = source->getHostMirror();
+    auto selectedMirror = selected.getHostMirror();
+    Kokkos::deep_copy(sourceMirror, source->getView());
+    Kokkos::deep_copy(selectedMirror, selected.getView());
+
+    const auto owned = source->getOwned();
+    const int nghost = source->getNghost();
+    nestedViewLoop(sourceMirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        std::array<ippl::Index, TestFixture::dim> pointIndices;
+        const std::array<size_t, TestFixture::dim> viewCoords{static_cast<size_t>(args)...};
+        for (unsigned d = 0; d < TestFixture::dim; ++d) {
+            const int global = owned[d].first()
+                               + static_cast<int>(viewCoords[d] - nghost) * owned[d].stride();
+            pointIndices[d] = ippl::Index(global, global);
+        }
+        auto point = std::make_from_tuple<ippl::NDIndex<TestFixture::dim>>(pointIndices);
+        const T value = sourceMirror(args...);
+        const T expected = subset.contains(point) && value > T(2) ? value : T(0);
+        assertEqual<T>(expected, selectedMirror(args...));
+    });
+}
+
+TYPED_TEST(BareFieldTest, SparseIndexedCompoundAssignment) {
+    using T = typename TestFixture::value_type;
+
+    auto& source = this->field;
+    auto view    = source->getView();
+    const auto lDom = source->getLayout().getLocalNDIndex();
+
+    Kokkos::parallel_for("Set field", source->getFieldRangePolicy(),
+                         FieldVal<TypeParam>{view, lDom});
+    Kokkos::fence();
+
+    typename TestFixture::field_type values(this->layout);
+    values = T(1);
+
+    ippl::SIndex<TestFixture::dim> sindex(this->layout);
+    sindex = *source > T(2);
+    values[sindex] = (*source)[sindex];
+    values[sindex] *= T(1.01);
+
+    auto sourceMirror = source->getHostMirror();
+    auto valuesMirror = values.getHostMirror();
+    Kokkos::deep_copy(sourceMirror, source->getView());
+    Kokkos::deep_copy(valuesMirror, values.getView());
+
+    const int nghost = source->getNghost();
+    nestedViewLoop(sourceMirror, nghost, [&]<typename... Idx>(const Idx... args) {
+        const T value = sourceMirror(args...);
+        const T expected = value > T(2) ? value * T(1.01) : T(1);
+        assertEqual<T>(expected, valuesMirror(args...));
+    });
+}
+
 TYPED_TEST(BareFieldTest, Sum) {
     using T = typename TestFixture::value_type;
 
diff --git a/unit_tests/CMakeLists.txt b/unit_tests/CMakeLists.txt
index af778d67b..90f555c0b 100644
--- a/unit_tests/CMakeLists.txt
+++ b/unit_tests/CMakeLists.txt
@@ -30,6 +30,7 @@ add_subdirectory(FEM)
 add_subdirectory(Interpolation)
 add_subdirectory(FFT)
 add_subdirectory(Field)
+add_subdirectory(Index)
 add_subdirectory(Meshes)
 add_subdirectory(Particle)
 add_subdirectory(PIC)
diff --git a/unit_tests/Field/Field.cpp b/unit_tests/Field/Field.cpp
index 505832b5d..f34f9c970 100644
--- a/unit_tests/Field/Field.cpp
+++ b/unit_tests/Field/Field.cpp
@@ -14,6 +14,9 @@
 template <typename>
 class FieldTest;
 
+template <typename>
+class Field3DRegressionTest;
+
 template <typename T, typename ExecSpace, unsigned Dim>
 class FieldTest<Parameters<T, ExecSpace, Rank<Dim>>> : public ::testing::Test {
 public:
@@ -63,6 +66,48 @@ class FieldTest<Parameters<T, ExecSpace, Rank<Dim>>> : public ::testing::Test {
     std::array<T, Dim> domain;
     };
 
+template <typename T, typename ExecSpace>
+class Field3DRegressionTest<Parameters<T, ExecSpace, Rank<3>>> : public ::testing::Test {
+public:
+    using value_type              = T;
+    using exec_space              = ExecSpace;
+    constexpr static unsigned dim = 3;
+    constexpr static size_t initial_size = 16;
+
+    using mesh_type      = ippl::UniformCartesian<T, dim>;
+    using centering_type = typename mesh_type::DefaultCentering;
+    using field_type     = ippl::Field<T, dim, mesh_type, centering_type, ExecSpace>;
+    using layout_type    = ippl::FieldLayout<dim>;
+
+    struct FieldBundle {
+        std::shared_ptr<layout_type> layout;
+        std::shared_ptr<mesh_type> mesh;
+        std::shared_ptr<field_type> field;
+    };
+
+    FieldBundle makeField(size_t n) const {
+        std::array<ippl::Index, dim> indices;
+        indices.fill(ippl::Index(n));
+        auto owned = std::make_from_tuple<ippl::NDIndex<dim>>(indices);
+
+        ippl::Vector<T, dim> hx;
+        ippl::Vector<T, dim> origin;
+        for (unsigned d = 0; d < dim; ++d) {
+            hx[d]     = T(1) / static_cast<T>(n);
+            origin[d] = T(0);
+        }
+
+        std::array<bool, dim> isParallel;
+        isParallel.fill(true);
+
+        FieldBundle bundle;
+        bundle.layout = std::make_shared<layout_type>(MPI_COMM_WORLD, owned, isParallel);
+        bundle.mesh   = std::make_shared<mesh_type>(owned, hx, origin);
+        bundle.field  = std::make_shared<field_type>(*bundle.mesh, *bundle.layout);
+        return bundle;
+    }
+};
+
 template <typename Params>
 struct VFieldVal {
     using vfield_view_type        = typename FieldTest<Params>::vfield_type::view_type;
@@ -151,6 +196,9 @@ struct FieldVal {
 using Tests = TestParams::tests<1, 2, 3, 4, 5, 6>;
 TYPED_TEST_SUITE(FieldTest, Tests);
 
+using Field3DRegressionTests = TestParams::tests<3>;
+TYPED_TEST_SUITE(Field3DRegressionTest, Field3DRegressionTests);
+
 TYPED_TEST(FieldTest, DeepCopy) {
     auto& field = this->field;
 
@@ -184,6 +232,33 @@ TYPED_TEST(FieldTest, Sum) {
     assertEqual<T>(expected, sum);
 }
 
+TYPED_TEST(Field3DRegressionTest, PowerOfTwoSubfieldSums) {
+    using T = typename TestFixture::value_type;
+    constexpr size_t N = TestFixture::initial_size;
+
+    auto source = this->makeField(N);
+    *source.field = T(1);
+
+    assertEqual<T>(static_cast<T>(N * N * N), source.field->sum());
+
+    for (size_t n = N / 2; n >= 2; n /= 2) {
+        auto target = this->makeField(N);
+        *target.field = T(0);
+
+        std::array<ippl::Index, TestFixture::dim> indices;
+        indices.fill(ippl::Index(n));
+        const auto subdomain = std::make_from_tuple<ippl::NDIndex<TestFixture::dim>>(indices);
+
+        (*target.field)[subdomain] = (*source.field)[subdomain];
+
+        assertEqual<T>(static_cast<T>(n * n * n), target.field->sum());
+
+        if (n == 2) {
+            break;
+        }
+    }
+}
+
 TYPED_TEST(FieldTest, Norm1) {
     using T = typename TestFixture::value_type;
 
diff --git a/unit_tests/Index/CMakeLists.txt b/unit_tests/Index/CMakeLists.txt
new file mode 100644
index 000000000..54cd13576
--- /dev/null
+++ b/unit_tests/Index/CMakeLists.txt
@@ -0,0 +1,4 @@
+file(RELATIVE_PATH _relPath "${PROJECT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}")
+message(STATUS "Adding unit tests found in ${_relPath}")
+
+add_ippl_test(Index USE_GTEST_MAIN)
diff --git a/unit_tests/Index/Index.cpp b/unit_tests/Index/Index.cpp
new file mode 100644
index 000000000..0fc0e11aa
--- /dev/null
+++ b/unit_tests/Index/Index.cpp
@@ -0,0 +1,46 @@
+//
+// Unit test Index
+//   Test strided Index semantics.
+//
+#include "Ippl.h"
+
+#include "gtest/gtest.h"
+
+TEST(IndexTest, StridedTouchesRequiresCommonElement) {
+    ippl::Index even(0, 10, 2);
+    ippl::Index odd(1, 9, 2);
+
+    EXPECT_FALSE(even.touches(odd));
+    EXPECT_FALSE(odd.touches(even));
+    EXPECT_TRUE(even.intersect(odd).empty());
+}
+
+TEST(IndexTest, StridedContainsRequiresContainedElements) {
+    ippl::Index even(0, 10, 2);
+    ippl::Index odd(1, 9, 2);
+    ippl::Index everyFourth(0, 8, 4);
+
+    EXPECT_FALSE(even.contains(odd));
+    EXPECT_TRUE(even.contains(everyFourth));
+    EXPECT_FALSE(everyFourth.contains(even));
+}
+
+TEST(IndexTest, GrowUsesIndexStrideForCoordinateBounds) {
+    ippl::Index even(2, 10, 2);
+    ippl::Index grown = even.grow(2);
+
+    EXPECT_EQ(grown.first(), -2);
+    EXPECT_EQ(grown.last(), 14);
+    EXPECT_EQ(grown.stride(), 2);
+    EXPECT_EQ(grown.length(), 9);
+}
+
+TEST(IndexTest, NegativeStrideGrowKeepsOrdering) {
+    ippl::Index descending(10, 2, -2);
+    ippl::Index grown = descending.grow(2);
+
+    EXPECT_EQ(grown.first(), 14);
+    EXPECT_EQ(grown.last(), -2);
+    EXPECT_EQ(grown.stride(), -2);
+    EXPECT_EQ(grown.length(), 9);
+}