diff --git a/src/FEM/GridPathSegmenter.h b/src/FEM/GridPathSegmenter.h
new file mode 100644
index 000000000..0e5f9a29a
--- /dev/null
+++ b/src/FEM/GridPathSegmenter.h
@@ -0,0 +1,29 @@
+#ifndef IPPL_GRIDPATH_SEGMENTER_H
+#define IPPL_GRIDPATH_SEGMENTER_H
+#include <array>
+
+namespace ippl {
+
+template<unsigned Dim, typename T>
+struct Segment {
+  ippl::Vector<T,Dim> p0, p1;
+};
+
+struct DefaultCellCrossingRule {};
+
+template<unsigned Dim, typename T, typename Rule = DefaultCellCrossingRule>
+struct GridPathSegmenter {
+
+  static KOKKOS_INLINE_FUNCTION
+  std::array<Segment<Dim,T>, Dim+1>
+  split(const ippl::Vector<T, Dim>& A,
+        const ippl::Vector<T, Dim>& B,
+        const ippl::Vector<T, Dim>& origin,
+        const ippl::Vector<T, Dim>& h);
+};
+
+} // namespace ippl 
+
+#include "GridPathSegmenter.hpp"
+
+#endif
diff --git a/src/FEM/GridPathSegmenter.hpp b/src/FEM/GridPathSegmenter.hpp
new file mode 100644
index 000000000..ea49aaf2f
--- /dev/null
+++ b/src/FEM/GridPathSegmenter.hpp
@@ -0,0 +1,123 @@
+#include <array>
+#include <cmath>
+#include <type_traits>
+
+namespace ippl {
+
+template<int I, int End, class F>
+KOKKOS_INLINE_FUNCTION
+constexpr void static_for(F&& f) {
+  if constexpr (I < End) {
+    f(std::integral_constant<int, I>{});
+    static_for<I+1, End>(std::forward<F>(f));
+  }
+}
+
+template<typename T>
+KOKKOS_INLINE_FUNCTION
+void sort2(T& a, T& b) { if (a > b) { T t=a; a=b; b=t; } }
+
+template<typename T>
+KOKKOS_INLINE_FUNCTION
+void sort3(T& a, T& b, T& c) { sort2(a,b); sort2(b,c); sort2(a,b); }
+
+template<unsigned Dim, typename T>
+KOKKOS_INLINE_FUNCTION
+ippl::Vector<T,Dim> lerp_point(const ippl::Vector<T,Dim>& A,
+                             const ippl::Vector<T,Dim>& B, T t) {
+  ippl::Vector<T,Dim> out{}; 
+  for (unsigned a=0; a<Dim; ++a) out[a] = A[a] + (B[a]-A[a]) * t;
+  return out;
+}
+
+// ---------- per-axis cut times (for DefaultCellCrossingRule) ----------
+template<unsigned Dim, typename T>
+struct CutTimes { std::array<T,Dim> t; };
+
+template<unsigned Dim, typename T>
+KOKKOS_INLINE_FUNCTION
+CutTimes<Dim,T> compute_axis_cuts_default(
+  const ippl::Vector<T,Dim>& A,
+  const ippl::Vector<T,Dim>& B,
+  const ippl::Vector<T,Dim>& origin,
+  const ippl::Vector<T,Dim>& h)
+{
+  CutTimes<Dim,T> cuts;
+  for (unsigned a=0; a<Dim; ++a) cuts.t[a] = (T)2; // sentinel (>1)
+
+  const T eps1 = (T)1e-12, one = (T)1;
+
+  auto axis_cut = [&](auto Ax) {
+    constexpr int a = Ax;
+    const T d   = B[a] - A[a];
+    const T eps = eps1 * h[a];
+    if (std::fabs(d) <= eps) return; // no motion → no crossing
+
+    const T nA  = (A[a] - origin[a]) / h[a]; // in cell units
+    // nearest plane index in direction toward B; bias off-plane a hair
+    const T k   = (d > 0) ? std::ceil(nA - eps1) : std::floor(nA + eps1);
+    const T pa  = origin[a] + k * h[a];
+    const T t   = (pa - A[a]) / d;
+
+    if (t > eps1 && t < one - eps1) cuts.t[a] = t;
+  };
+
+  static_for<0,Dim>(axis_cut);
+  return cuts;
+}
+
+// ---------- endpoints builder: [0, t_sorted..., 1] (keeps duplicates) ----------
+template<unsigned Dim, typename T>
+KOKKOS_INLINE_FUNCTION
+void make_endpoints_fixed(const CutTimes<Dim,T>& cuts,
+                          std::array<T,Dim+2>& Tcuts /*out*/)
+{
+  T t0 = cuts.t[0];
+  T t1 = (Dim>=2) ? cuts.t[1] : (T)2;
+  T t2 = (Dim==3) ? cuts.t[2] : (T)2;
+
+  if constexpr (Dim==2) {
+    sort2(t0, t1);
+  } else { // Dim==3
+    sort3(t0, t1, t2);
+  }
+
+  const T one = (T)1;
+  auto clamp_or_one = [&](T v)->T { return (v >= (T)1.5) ? one : (v < one ? v : one); };
+
+  Tcuts[0] = (T)0;
+  Tcuts[1] = clamp_or_one(t0);
+  if constexpr (Dim>=2) Tcuts[2] = clamp_or_one(t1);
+  if constexpr (Dim==3) Tcuts[3] = clamp_or_one(t2);
+  Tcuts[Dim+1] = one;
+}
+
+// ---------------------------------------------------------------------------------
+// DefaultCellCrossingRule specialization 
+// ---------------------------------------------------------------------------------
+
+template<unsigned Dim, typename T, typename Rule>
+KOKKOS_INLINE_FUNCTION
+std::array<Segment<Dim,T>, Dim+1>
+GridPathSegmenter<Dim,T,Rule>::split(
+    const ippl::Vector<T,Dim>& A,
+    const ippl::Vector<T,Dim>& B,
+    const ippl::Vector<T,Dim>& origin,
+    const ippl::Vector<T,Dim>& h)
+{
+  const auto cuts = compute_axis_cuts_default<Dim,T>(A,B,origin,h);
+
+  std::array<T,Dim+2> Tcuts{};
+  make_endpoints_fixed<Dim,T>(cuts, Tcuts);
+
+  std::array<Segment<Dim,T>, Dim+1> segs{};
+  for (unsigned i = 0; i < Dim + 1; ++i) {
+    const T ta = Tcuts[i];
+    const T tb = Tcuts[i+1];
+    segs[i].p0 = lerp_point<Dim,T>(A,B,ta);
+    segs[i].p1 = lerp_point<Dim,T>(A,B,tb);
+  }
+  return segs;
+}
+
+} // namespace ippl 
diff --git a/src/FEM/ProjectCurrent.hpp b/src/FEM/ProjectCurrent.hpp
new file mode 100644
index 000000000..b4970fbcf
--- /dev/null
+++ b/src/FEM/ProjectCurrent.hpp
@@ -0,0 +1,42 @@
+#ifndef IPPL_PROJECT_CURRENT_H
+#define IPPL_PROJECT_CURRENT_H
+
+namespace ippl {
+
+
+template <typename Mesh,
+          typename ChargeAttrib,
+          typename PosAttrib,
+          typename FEMVector,
+          typename NedelecSpace,
+          typename policy_type = Kokkos::RangePolicy<>>
+inline void assemble_current_whitney1(const Mesh& mesh,
+                                      const ChargeAttrib& q_attrib,
+                                      const PosAttrib& X0,
+                                      const PosAttrib& X1,
+                                      FEMVector& fem_vector,
+                                      const NedelecSpace& space,
+                                      policy_type iteration_policy)
+{
+  using T = Mesh::value_type;
+
+  const auto origin = mesh.getOrigin();
+  const auto h   = mesh.getMeshSpacing();
+  constexpr unsigned Dim = Mesh::Dimension;
+
+  Kokkos::parallel_for("assemble_current_whitney1_make_segments", iteration_policy,
+    KOKKOS_LAMBDA(const std::size_t p) {
+
+    auto segs = ippl::GridPathSegmenter<Dim, T, ippl::DefaultCellCrossingRule>
+                    ::split(X0(p), X1(p), origin, h);
+    
+     
+
+    
+  });
+
+
+}
+
+}
+#endif
diff --git a/src/IpplCore.h b/src/IpplCore.h
index 4e21f50c6..4a1b750f2 100644
--- a/src/IpplCore.h
+++ b/src/IpplCore.h
@@ -56,5 +56,7 @@
 
 // FEM Operators
 #include "FEM/FEMInterpolate.hpp"
+#include "FEM/GridPathSegmenter.h"
+#include "FEM/ProjectCurrent.hpp"
 
 #endif
diff --git a/unit_tests/FEM/AssembleCurrentRHS.cpp b/unit_tests/FEM/AssembleCurrentRHS.cpp
new file mode 100644
index 000000000..dfdd6ef03
--- /dev/null
+++ b/unit_tests/FEM/AssembleCurrentRHS.cpp
@@ -0,0 +1,171 @@
+#include "Ippl.h"
+#include "TestUtils.h"
+#include "gtest/gtest.h"
+
+template <class PLayout>
+struct Bunch : public ippl::ParticleBase<PLayout> {
+    Bunch(PLayout& playout)
+        : ippl::ParticleBase<PLayout>(playout) {
+        this->addAttribute(Q);
+    }
+
+    ~Bunch() {}
+
+    typedef ippl::ParticleAttrib<double> charge_container_type;
+    charge_container_type Q;
+};
+
+
+
+template <unsigned Dim, typename T> struct ElemSelector;
+
+template <typename T>
+struct ElemSelector<1, T> {
+  using Elem = ippl::EdgeElement<T>;
+  using Quad = ippl::MidpointQuadrature<T, 1, Elem>;
+
+  static Elem make_elem() { return Elem{}; }
+  static Quad make_quad(const Elem& e) { return Quad(e); }
+};
+
+template <typename T>
+struct ElemSelector<2, T> {
+  using Elem = ippl::QuadrilateralElement<T>;
+  using Quad = ippl::MidpointQuadrature<T, 1, Elem>;
+
+  static Elem make_elem() { return Elem{}; }
+  static Quad make_quad(const Elem& e) { return Quad(e); }
+};
+
+template <typename T>
+struct ElemSelector<3, T> {
+  using Elem = ippl::HexahedralElement<T>;
+  using Quad = ippl::MidpointQuadrature<T, 1, Elem>;
+
+  static Elem make_elem() { return Elem{}; }
+  static Quad make_quad(const Elem& e) { return Quad(e); }
+};
+
+
+template <typename>
+class AssembleCurrentTest;
+
+template <typename T, unsigned Dim>
+class AssembleCurrentTest<Parameters<T, Rank<Dim>>> : public ::testing::Test {
+public:
+  using value_type = T;
+  static constexpr unsigned dim = Dim;
+
+  using Mesh_t      = ippl::UniformCartesian<T, dim>;
+  using Centering_t = typename Mesh_t::DefaultCentering;
+
+  using ElemSel     = ElemSelector<dim, T>;
+  using Elem        = typename ElemSel::Elem;
+  using Quad        = typename ElemSel::Quad;
+  // Somehow the field type is completely irrelevant for the fem space??
+  using FieldType   = ippl::Field<T, Dim, Mesh_t, typename Mesh_t::DefaultCentering>;
+  using Layout      = ippl::FieldLayout<Dim>;
+
+  using NedelecType = ippl::NedelecSpace<T, Dim, 1, Elem, Quad, FieldType>;
+
+  using playout_t   = ippl::ParticleSpatialLayout<T, dim>;
+  using bunch_t     = Bunch<playout_t>;
+
+  static ippl::NDIndex<dim> make_owned_nd(int nx) {
+    ippl::Index I0(nx);
+    if constexpr (dim == 1)
+      return ippl::NDIndex<1>(I0);
+    else if constexpr (dim == 2)
+      return ippl::NDIndex<2>(I0, I0);
+    else
+      return ippl::NDIndex<3>(I0, I0, I0);
+  }
+
+  static Layout make_layout(const ippl::NDIndex<dim>& owned) {
+    std::array<bool, dim> par{}; par.fill(true);
+    return ippl::FieldLayout<dim>(MPI_COMM_WORLD, owned, par);
+  }
+
+  static Mesh_t make_mesh(const ippl::NDIndex<dim>& owned,
+                          const ippl::Vector<T, dim>& h,
+                          const ippl::Vector<T, dim>& origin) {
+    return Mesh_t(owned, h, origin);
+  }
+
+  static NedelecType make_space(Mesh_t& mesh, Layout l) {
+    Elem e = ElemSel::make_elem();
+    Quad q = ElemSel::make_quad(e);
+    return NedelecType(mesh, e, q, l);
+  }
+};
+
+using Precisions = TestParams::Precisions;
+using Ranks      = TestParams::Ranks<2, 3>;
+using Tests      = TestForTypes<CreateCombinations<Precisions, Ranks>::type>::type;
+
+TYPED_TEST_SUITE(AssembleCurrentTest, Tests);
+
+// ------------------ Actual minimal test ------------------
+
+TYPED_TEST(AssembleCurrentTest, SingleParticle_Smoke) {
+  using T   = typename TestFixture::value_type;
+  constexpr unsigned Dim = TestFixture::dim;
+
+  using bunch_t    = typename TestFixture::bunch_t;
+  using playout_t  = typename TestFixture::playout_t;
+
+  int nx = 4;
+
+  ippl::Vector<T, Dim> origin(0.0);
+  ippl::Vector<T, Dim> h(1.0);
+
+  auto owned  = TestFixture::make_owned_nd(nx);
+  auto layout = TestFixture::make_layout(owned);
+  auto mesh   = TestFixture::make_mesh(owned, h, origin);
+  auto space  = TestFixture::make_space(mesh, layout);
+
+  playout_t playout(layout, mesh);
+  bunch_t   bunch(playout);
+
+  // --- create 1 particle ---
+  bunch.create(1);
+  {
+    auto R_host = bunch.R.getHostMirror();
+    auto Q_host = bunch.Q.getHostMirror();
+    for (unsigned d=0; d<Dim; ++d)
+      R_host(0)[d] = T(0.5) * h[d];  // inside domain
+    Q_host(0) = 1.0;
+    Kokkos::deep_copy(bunch.R.getView(), R_host);
+    Kokkos::deep_copy(bunch.Q.getView(), Q_host);
+    bunch.update();
+  }
+
+  // --- mock "next" positions (slightly moved) ---
+  bunch_t bunch_next(playout);
+  bunch_next.create(1);
+  {
+    auto Rn_host = bunch_next.R.getHostMirror();
+    for (unsigned d=0; d<Dim; ++d)
+      Rn_host(0)[d] = T(0.6) * h[d];
+    Kokkos::deep_copy(bunch_next.R.getView(), Rn_host);
+    bunch_next.update();
+  }
+
+  auto policy = Kokkos::RangePolicy<>(0, bunch.getLocalNum());
+  auto fem_vector = space.createFEMVector();
+
+  // Call function (no assertions yet)
+  ippl::assemble_current_whitney1(mesh, bunch.Q, bunch.R, bunch_next.R, fem_vector, space, policy);
+
+
+  SUCCEED() << "assemble_current_whitney1 ran without error for 1 particle in "
+            << Dim << "D.";
+}
+
+int main(int argc, char* argv[]) {
+  ippl::initialize(argc, argv);
+  ::testing::InitGoogleTest(&argc, argv);
+  int result = RUN_ALL_TESTS();
+  ippl::finalize();
+  return result;
+}
diff --git a/unit_tests/FEM/CMakeLists.txt b/unit_tests/FEM/CMakeLists.txt
index 61f614e9b..b15a1273f 100644
--- a/unit_tests/FEM/CMakeLists.txt
+++ b/unit_tests/FEM/CMakeLists.txt
@@ -1,6 +1,8 @@
 file(RELATIVE_PATH _relPath "${PROJECT_SOURCE_DIR}" "${CMAKE_CURRENT_SOURCE_DIR}")
 message(STATUS "Adding unit tests found in ${_relPath}")
 
+add_ippl_test(AssembleCurrentRHS)
+add_ippl_test(GenerateSegments)
 add_ippl_test(EdgeElement)
 add_ippl_test(QuadrilateralElement)
 add_ippl_test(HexahedralElement)
diff --git a/unit_tests/FEM/GenerateSegments.cpp b/unit_tests/FEM/GenerateSegments.cpp
new file mode 100644
index 000000000..fb6ce4526
--- /dev/null
+++ b/unit_tests/FEM/GenerateSegments.cpp
@@ -0,0 +1,235 @@
+#include <iostream>
+#include <gtest/gtest.h>
+
+#include "Ippl.h"
+#include "TestUtils.h"
+
+
+template <typename>
+class GridPathSegmenterTest;
+
+template <typename T, unsigned Dim>
+class GridPathSegmenterTest<Parameters<T, Rank<Dim>>> : public ::testing::Test {
+public:
+  using value_type = T;
+  static constexpr unsigned dim = Dim;
+
+  using Rule = ippl::DefaultCellCrossingRule;
+
+  struct Scenario {
+    ippl::Vector<T,Dim> A{};
+    ippl::Vector<T,Dim> B{};
+    std::size_t expected_segments = 0; // expected REAL segments (after compacting zero-length)
+    const char* name = "";
+  };
+
+private:
+
+  static KOKKOS_INLINE_FUNCTION T plane(unsigned a, std::size_t k,
+                                        const ippl::Vector<T,Dim>& origin,
+                                        const ippl::Vector<T,Dim>& h) {
+    return origin[a] + T(k) * h[a];
+  }
+
+  static ippl::Vector<T,Dim> basis_scaled(unsigned a, const ippl::Vector<T,Dim>& h, T s) {
+    ippl::Vector<T,Dim> v{}; 
+    for (unsigned d=0; d<Dim; ++d) v[d]=T(0);
+    v[a] = s * h[a];
+    return v;
+  }
+
+public:
+
+  // Same cell: no plane crossings → 1 real segment
+  Scenario same_cell(const ippl::Vector<T,Dim>& origin,
+                     const ippl::Vector<T,Dim>& h) const {
+    ippl::Vector<T,Dim> A{}, B{};
+    for (unsigned d=0; d<Dim; ++d) {
+      A[d] = origin[d] + T(0.25) * h[d];
+      B[d] = origin[d] + T(0.75) * h[d];
+    }
+    return {A, B, 1u, "same_cell"};
+  }
+
+  // Single-axis forward cut: cross plane once along +axis → 2 segments
+  Scenario single_axis_forward(unsigned axis,
+                               const ippl::Vector<T,Dim>& origin,
+                               const ippl::Vector<T,Dim>& h) const {
+    ippl::Vector<T,Dim> A{}, B{};
+    for (unsigned d=0; d<Dim; ++d) A[d] = origin[d] + T(0.2) * h[d];
+    // Start near the upper face along 'axis' to guarantee a crossing
+    A[axis] = origin[axis] + T(0.9) * h[axis];
+    B = A;
+    B[axis] += T(0.25) * h[axis]; // crosses plane at k+1 once
+    return {A, B, 2u, "single_axis_forward"};
+  }
+
+  Scenario two_axes_forward(unsigned axis0, unsigned axis1,
+                            const ippl::Vector<T,Dim>& origin,
+                            const ippl::Vector<T,Dim>& h) const {
+    static_assert(Dim >= 2, "two_axes_forward requires Dim>=2");
+    ippl::Vector<T,Dim> A{}, B{};
+    for (unsigned d=0; d<Dim; ++d) A[d] = origin[d] + T(0.2) * h[d];
+    A[axis0] = origin[axis0] + T(0.9) * h[axis0];
+    A[axis1] = origin[axis1] + T(0.85) * h[axis1];
+    B = A;
+    B[axis0] += T(0.25) * h[axis0];
+    B[axis1] += T(0.30) * h[axis1];
+    return {A, B, 3u, "two_axes_forward"};
+  }
+
+  Scenario three_axes_forward(const ippl::Vector<T,Dim>& origin,
+                              const ippl::Vector<T,Dim>& h) const {
+    static_assert(Dim >= 3, "three_axes_forward requires Dim>=3");
+    ippl::Vector<T,Dim> A{}, B{};
+
+    A[0] = origin[0] + T(0.9)  * h[0];  // dist to plane = 0.1h
+    A[1] = origin[1] + T(0.85) * h[1];  // dist to plane = 0.15h
+    A[2] = origin[2] + T(0.8)  * h[2];  // dist to plane = 0.20h
+
+    B = A;
+
+    for (unsigned d = 0; d < 3; ++d) {
+      B[d] += T(0.30) * h[d];
+    }
+    return {A, B, 4u, "three_axes_forward"};
+  }
+
+  Scenario start_on_plane(unsigned axis,
+                          const ippl::Vector<T,Dim>& origin,
+                          const ippl::Vector<T,Dim>& h) const {
+    ippl::Vector<T,Dim> A{}, B{};
+    for (unsigned d=0; d<Dim; ++d) A[d] = origin[d] + T(0.3) * h[d];
+    // Put A exactly at plane k=1 on 'axis'
+    A[axis] = plane(axis, 1, origin, h);
+    B = A;
+    B[axis] += T(0.2) * h[axis];
+    return {A, B, 1u, "start_on_plane"};
+  }
+
+  Scenario vertex_hit(const ippl::Vector<T,Dim>& origin,
+                      const ippl::Vector<T,Dim>& h) const {
+    ippl::Vector<T,Dim> A{}, B{};
+    for (unsigned d = 0; d < Dim; ++d) {
+      A[d] = origin[d] + T(0.9) * h[d];
+    }
+    B = A;
+    for (unsigned d = 0; d < Dim; ++d) {
+      B[d] += T(0.2) * h[d];
+    }
+    return {A, B, 2u, "vertex_hit"};
+  }
+
+  static ippl::Vector<T,Dim> ones() {
+    ippl::Vector<T,Dim> a{}; for (unsigned d=0; d<Dim; ++d) a[d] = T(1); return a;
+  }
+  static ippl::Vector<T,Dim> zeros() {
+    ippl::Vector<T,Dim> a{}; for (unsigned d=0; d<Dim; ++d) a[d] = T(0); return a;
+  }
+
+  struct OriginH {
+    ippl::Vector<T,Dim> origin{};
+    ippl::Vector<T,Dim> h{};
+    const char* name{};
+  };
+
+  std::vector<Scenario>
+  scenario_cases(const ippl::Vector<T,Dim>& origin, const ippl::Vector<T,Dim>& h) {
+    std::vector<Scenario> v;
+
+    // Always useful
+    v.push_back(same_cell(origin, h));
+    v.push_back(vertex_hit(origin, h));
+    v.push_back(single_axis_forward(0, origin, h));
+    v.push_back(start_on_plane(0, origin, h));
+
+
+    // Dim >= 2 scenarios
+    if constexpr (Dim >= 2) {
+      v.push_back(two_axes_forward(0, 1, origin, h));
+    }
+
+    // Dim >= 3 scenarios
+    if constexpr (Dim >= 3) {
+      v.push_back(three_axes_forward(origin, h));
+    }
+
+    return v;
+  }
+
+  static std::vector<OriginH> origin_h_cases() {
+    return std::vector{ // C++20 CTAD for vectors
+      OriginH{ /*origin=*/fill_val(T(0)), /*h=*/fill_val(T(1)), "origin=0, h=1" },
+      OriginH{ /*origin=*/fill_seq(T(-0.7), T(0.3)), /*h=*/fill_val(T(1)), "shifted origin, h=1" },
+      OriginH{ /*origin=*/fill_val(T(0)), /*h=*/fill_seq(T(0.3), T(0.7)),"origin=0, nonunit h" },
+      OriginH{ /*origin=*/fill_seq(T(-1.1), T(0.5)), /*h=*/fill_seq(T(0.25), T(1.1)), "shifted origin, mixed h" },
+    };
+  }
+
+private:
+  static ippl::Vector<T,Dim> fill_val(T v) {
+    ippl::Vector<T,Dim> a; 
+    a = v; 
+    return a;
+  }
+  static ippl::Vector<T,Dim> fill_seq(T base, T step) {
+    ippl::Vector<T,Dim> a{};
+    for (unsigned d=0; d<Dim; ++d) a[d] = base + step*T(d);
+    return a;
+  }
+};
+
+using Precisions = TestParams::Precisions;
+using Ranks      = TestParams::Ranks<1, 2, 3>;
+using Tests      = TestForTypes<CreateCombinations<Precisions, Ranks>::type>::type;
+
+TYPED_TEST_SUITE(GridPathSegmenterTest, Tests);
+
+TYPED_TEST(GridPathSegmenterTest, SegmentCountMatchesScenario) {
+  using T   = typename TestFixture::value_type;
+  constexpr unsigned Dim = TestFixture::dim;
+
+  auto seg_is_real = [](const auto& s) {
+    T maxc = T(0);
+    for (unsigned d=0; d<Dim; ++d) {
+      maxc = std::max({maxc, std::abs(s.p0[d]), std::abs(s.p1[d])});
+    }
+    const T eps = std::numeric_limits<T>::epsilon() * (T)100; // a bit looser than ulp
+    for (unsigned d=0; d<Dim; ++d) {
+      if (std::abs(s.p1[d] - s.p0[d]) > eps * (T(1) + maxc)) return true;
+    }
+    return false;
+  };
+
+  for (const auto& oh : TestFixture::origin_h_cases()) {
+    auto scenarios = TestFixture::scenario_cases(oh.origin, oh.h);
+    for (const auto& sc : scenarios) {
+      // run split
+      auto segs = ippl::GridPathSegmenter<Dim, T, typename TestFixture::Rule>
+                    ::split(sc.A, sc.B, oh.origin, oh.h);
+
+      // compact: count only real (non-zero-length) segments
+      std::size_t real_count = 0;
+      for (std::size_t i = 0; i < Dim + 1; ++i) {
+        if (seg_is_real(segs[i])) ++real_count;
+      }
+
+      EXPECT_EQ(real_count, sc.expected_segments)
+        << "Mismatch in segment count for scenario='" << sc.name
+        << "' with origin/h case='" << oh.name << "'\n"
+        << "  expected=" << sc.expected_segments
+        << "  got=" << real_count;
+    }
+  }
+}
+
+int main(int argc, char* argv[]) {
+  ippl::initialize(argc, argv);
+  int result = 0;
+  {
+    ::testing::InitGoogleTest(&argc, argv);
+    result = RUN_ALL_TESTS();
+  }
+  ippl::finalize();
+  return result;
+}