Factor out a Neighbors to package up all 6 neighbors when updating cells

This will allow more flexibility in boundary conditions.

It seems to have similar perf as passing only the 3 cells within about 2%.

src/sim.rs

@@ -98,10 +98,7 @@ impl<M: Material + Clone + Default + Send + Sync> SimState<M> {
         Zip::from(ndarray::indices_of(&self.cells)).and(&self.cells).par_apply_assign_into(
             &mut scratch_cells,
             |(z, y, x), cell| {
-                let down_cell = self.cells.get([z, y+1, x]);
+                cell.clone().step_h(Neighbors::new(&self.cells, [z, y, x]), half_time_step.into(), feature_size.into())
-                let right_cell = self.cells.get([z, y, x+1]);
-                let in_cell = self.cells.get([z+1, y, x]);
-                cell.clone().step_h(right_cell, down_cell, in_cell, half_time_step.into(), feature_size.into())
             });
         trace!("step_h end");
 
@@ -110,10 +107,7 @@ impl<M: Material + Clone + Default + Send + Sync> SimState<M> {
         Zip::from(ndarray::indices_of(&scratch_cells)).and(&scratch_cells).par_apply_assign_into(
             &mut self.cells,
             |(z, y, x), cell| {
-                let left_cell = scratch_cells.get([z, y, x.wrapping_sub(1)]);
+                cell.clone().step_e(Neighbors::new(&scratch_cells, [z, y, x]), half_time_step.into(), feature_size.into())
-                let up_cell = scratch_cells.get([z, y.wrapping_sub(1), x]);
-                let out_cell = scratch_cells.get([z.wrapping_sub(1), y, x]);
-                cell.clone().step_e(left_cell, up_cell, out_cell, half_time_step.into(), feature_size.into())
             });
         trace!("step_e end");
 
@@ -247,6 +241,30 @@ pub struct Cell<M = mat::Static> {
     mat: M,
 }
 
+struct Neighbors<'a, Cell> {
+    left: Option<&'a Cell>,
+    right: Option<&'a Cell>,
+    up: Option<&'a Cell>,
+    down: Option<&'a Cell>,
+    out: Option<&'a Cell>,
+    in_: Option<&'a Cell>,
+}
+
+impl<'a, Cell> Neighbors<'a, Cell> {
+    fn new(array: &'a Array3<Cell>, location: [usize; 3]) -> Self {
+        let [z, y, x] = location;
+        let left  = array.get([z, y, x.wrapping_sub(1)]);
+        let right = array.get([z, y, x + 1]);
+        let down = array.get([z, y.wrapping_sub(1), x]);
+        let up = array.get([z, y + 1, x]);
+        let out = array.get([z.wrapping_sub(1), y, x]);
+        let in_ = array.get([z + 1, y, x]);
+        Self {
+            left, right, up, down, out, in_
+        }
+    }
+}
+
 impl<M> Cell<M> {
     pub fn ex(&self) -> Flt {
         self.state.e().x()
@@ -297,7 +315,7 @@ impl<M: Material> Cell<M> {
         self.state.h = self.mat.step_h(&self.state, Vec3::new(0.0, 0.0, delta_bz));
     }
 
-    fn step_h(mut self, right: Option<&Self>, down: Option<&Self>, in_: Option<&Self>, delta_t: Real, feature_size: Real) -> Self {
+    fn step_h(mut self, neighbors: Neighbors<'_, Self>, delta_t: Real, feature_size: Real) -> Self {
         // ```tex
         // Maxwell-Faraday equation: $\nabla x E = -dB/dt$
         //
@@ -313,21 +331,21 @@ impl<M: Material> Cell<M> {
         use consts::real::ZERO;
         let inv_feature_size = Real::from_inner(1.0) / feature_size;
 
-        let (delta_ey_x, delta_ez_x) = match right {
+        let (delta_ey_x, delta_ez_x) = match neighbors.right {
             None => (ZERO(), ZERO()),
             Some(right) => (
                 inv_feature_size * (right.ey() - self.ey()),
                 inv_feature_size * (right.ez() - self.ez()),
             )
         };
-        let (delta_ex_y, delta_ez_y) = match down {
+        let (delta_ex_y, delta_ez_y) = match neighbors.down {
             None => (ZERO(), ZERO()),
             Some(down) => (
                 inv_feature_size * (down.ex() - self.ex()),
                 inv_feature_size * (down.ez() - self.ez()),
             )
         };
-        let (delta_ex_z, delta_ey_z) = match in_ {
+        let (delta_ex_z, delta_ey_z) = match neighbors.in_ {
             None => (ZERO(), ZERO()),
             Some(in_) => (
                 inv_feature_size * (in_.ex() - self.ex()),
@@ -349,7 +367,7 @@ impl<M: Material> Cell<M> {
 
     /// delta_t = timestep covered by this function. i.e. it should be half the timestep of the simulation
     /// feature_size = how many units apart is the center of each adjacent cell on the grid.
-    fn step_e(self, left: Option<&Self>, up: Option<&Self>, out: Option<&Self>, delta_t: Real, feature_size: Real) -> Self {
+    fn step_e(self, neighbors: Neighbors<'_, Self>, delta_t: Real, feature_size: Real) -> Self {
         // ```tex
         // Ampere's circuital law with Maxwell's addition, in SI units ("macroscopic version"):
         //   $\nabla x H = J_f + dD/dt$    where $J_f$ = current density = $\sigma E$, $\sigma$ being a material parameter ("conductivity")
@@ -371,6 +389,12 @@ impl<M: Material> Cell<M> {
         // | d/dx d/dy d/dz | = | dHx/dz - dHz/dx |
         // |                |   |                 |
         // |  Hx   Hy   Hz  |   [ dHy/dx - dHx/dy ]
+        //
+        // Boundary conditions (TODO):
+        // For something on the yz plane, E(x=0, t+1) should be E(x=1, t). i.e. the wave travels
+        // delta_x = 1 per timestep. Roughly, we want to take the line from E(x=0, t), E(x=1, t)
+        // and continue that to E(x=-1, t) ?
+        //
         // ```
         
         use consts::real::{EPS0, TWO, ZERO};
@@ -378,21 +402,21 @@ impl<M: Material> Cell<M> {
         let sigma: Vec3 = self.mat.conductivity();
         let inv_feature_size = Real::from_inner(1.0) / feature_size;
 
-        let (delta_hy_x, delta_hz_x) = match left {
+        let (delta_hy_x, delta_hz_x) = match neighbors.left {
             None => (ZERO(), ZERO()),
             Some(left) => (
                 inv_feature_size * (self.hy() - left.hy()),
                 inv_feature_size * (self.hz() - left.hz()),
             )
         };
-        let (delta_hx_y, delta_hz_y) = match up {
+        let (delta_hx_y, delta_hz_y) = match neighbors.up {
             None => (ZERO(), ZERO()),
             Some(up) => (
                 inv_feature_size * (self.hx() - up.hx()),
                 inv_feature_size * (self.hz() - up.hz()),
             )
         };
-        let (delta_hx_z, delta_hy_z) = match out {
+        let (delta_hx_z, delta_hy_z) = match neighbors.out {
             None => (ZERO(), ZERO()),
             Some(out) => (
                 inv_feature_size * (self.hx() - out.hx()),