new purely functional iterator

include/Eisvogel/Interpolator.hh

@@ -3,12 +3,50 @@
 
 #include <cmath>
 #include <iostream>
+#include <type_traits>
 
 #include "Common.hh"
 #include "Kernels.hh"
 #include "NDArray.hh"
 #include "IteratorUtils.hh"
 
+#include "CoordUtils.hh"
+
+template <typename FuncT, typename KernelT,
+	  typename ValueT = std::invoke_result_t<FuncT, GridVector&>>
+ValueT InterpolateFunc(FuncT func, KernelT& kernel, CoordVector& target_inds) {
+
+  std::size_t dims = target_inds.size();
+
+  GridVector start_inds(dims, 0);
+  GridVector end_inds(dims, 0);
+
+  for(std::size_t i = 0; i < dims; i++) {
+    scalar_t start_ind = std::ceil(target_inds(i) - kernel.Support());
+    scalar_t end_ind = std::floor(target_inds(i) + kernel.Support() + 1);
+
+    start_inds(i) = int(start_ind);
+    end_inds(i) = int(end_ind);
+  }
+
+  ValueT interpolated_value = ValueT();
+
+  // iterate over all dimensions
+  for(GridCounter cnt(start_inds, end_inds); cnt.running(); ++cnt) {
+
+    scalar_t kernel_weight = 1.0;
+    for(std::size_t i = 0; i < dims; i++) {
+      kernel_weight *= kernel(target_inds(i) - cnt(i));
+    }
+
+    ValueT cur_val = func(cnt.index());
+
+    interpolated_value += cur_val * kernel_weight;      
+  }
+
+  return interpolated_value;
+}
+
 template <template<typename, std::size_t> class ArrayT, typename ValueT, std::size_t dims>
 class Interpolator {
 

include/Eisvogel/IteratorUtils.hh

@@ -40,6 +40,44 @@ public:
   inline IndexVector& index() {return m_cur;}
 };
 
+
+class GridCounter {
+
+private:
+
+  GridVector m_start, m_end, m_cur;
+  std::size_t number_dims;
+
+public:
+
+  GridCounter(const GridVector& start, const GridVector& end) : m_start(start), m_end(end), m_cur(start), 
+								   number_dims(start.size()) { }
+  GridCounter& operator++() {
+    for(size_t i = 0; i < number_dims; i++) {
+      if(++m_cur(i) == m_end(i)) {
+	if(i < number_dims - 1) {
+	  m_cur(i) = m_start(i);
+	}
+      }
+      else
+	break;
+    }
+
+    return *this;
+  }
+
+  bool running() {
+    return m_cur(number_dims - 1) < m_end(number_dims - 1);
+  }
+
+  GridVector::type& operator()(size_t ind) {return m_cur(ind);}
+
+  auto begin() {return m_cur.begin();}
+  auto end() {return m_cur.end();}
+
+  inline GridVector& index() {return m_cur;}
+};
+
 // some additional utility functions
 bool isInIndexRange(const IndexVector& inds, const IndexVector& start_inds, const IndexVector& stop_inds);
 

include/Eisvogel/NDArray.hh

@@ -9,6 +9,8 @@
 
 #include "Serialization.hh"
 
+// TODO: probably also want a fixed-sized version of NDArray that allows to specify the individual dimensions
+
 template <class T, std::size_t dims>
 class NDArray {
 
@@ -277,6 +279,7 @@ template <class T>
 using DenseVector = DenseNDArray<T, 1>;
 
 using IndexVector = DenseVector<std::size_t>;
+using GridVector = DenseVector<unsigned int>;
 
 template <class T>
 using ScalarField3D = DenseNDArray<T, 3>;

src/Integrator.cxx

@@ -65,10 +65,64 @@ scalar_t Integrator::integrate(scalar_t t, const Current0D& curr, scalar_t os_fa
 
       CoordVector wf_eval_frac_inds = m_dwf -> getFracInds(wf_eval_pos);
 
-      FieldVector wf_rzphi = CU::MakeFieldVectorRZPHI(m_itpl_E_r -> Interpolate(wf_eval_frac_inds),
-						      m_itpl_E_z -> Interpolate(wf_eval_frac_inds),
-						      m_itpl_E_phi -> Interpolate(wf_eval_frac_inds));
-
+      auto get_E_r = [&](GridVector& vec) -> scalar_t {
+
+	int ind_t = vec(0);
+	int ind_z = vec(1);
+	int ind_r = vec(2);	
+
+	if(ind_t < 0) {
+	  return 0.0;
+	}
+
+	if(ind_r < 0) {
+	  ind_r = -ind_r;
+	}
+
+	IndexVector ind = {(std::size_t)ind_t, (std::size_t)ind_z, (std::size_t)ind_r};	
+	return m_dwf -> E_r().operator()(ind);
+      };
+
+      auto get_E_z = [&](GridVector& vec) -> scalar_t {
+
+	int ind_t = vec(0);
+	int ind_z = vec(1);
+	int ind_r = vec(2);	
+
+	if(ind_t < 0) {
+	  return 0.0;
+	}
+
+	if(ind_r < 0) {
+	  ind_r = -ind_r;
+	}
+
+	IndexVector ind = {(std::size_t)ind_t, (std::size_t)ind_z, (std::size_t)ind_r};	
+	return m_dwf -> E_z().operator()(ind);
+      };
+
+      auto get_E_phi = [&](GridVector& vec) -> scalar_t {
+
+	int ind_t = vec(0);
+	int ind_z = vec(1);
+	int ind_r = vec(2);	
+
+	if(ind_t < 0) {
+	  return 0.0;
+	}
+
+	if(ind_r < 0) {
+	  ind_r = -ind_r;
+	}
+
+	IndexVector ind = {(std::size_t)ind_t, (std::size_t)ind_z, (std::size_t)ind_r};	
+	return m_dwf -> E_phi().operator()(ind);
+      };
+
+      FieldVector wf_rzphi = CU::MakeFieldVectorRZPHI(InterpolateFunc(get_E_r, *m_kernel, wf_eval_frac_inds),
+						      InterpolateFunc(get_E_z, *m_kernel, wf_eval_frac_inds),
+						      InterpolateFunc(get_E_phi, *m_kernel, wf_eval_frac_inds));
+
       FieldVector wf_xyz = CU::RZPHI_to_XYZ(wf_rzphi, cur_pos_txyz);
 
       scalar_t wf_val = CU::getXComponent(wf_xyz) * CU::getX(segment_velocity) +