From 0b0d03541df161177a92a27c05eadd01aeba166b Mon Sep 17 00:00:00 2001
From: aprilnovak <novak@berkeley.edu>
Date: Thu, 15 Oct 2020 16:16:34 -0500
Subject: [PATCH 1/8] Put get_indices implementation into StructuredMesh base
 class. Refs #1695

---
 include/openmc/mesh.h | 21 ++++++++++++---------
 src/mesh.cpp          | 35 +++++++++++++++--------------------
 2 files changed, 27 insertions(+), 29 deletions(-)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index b7dcb017b4..75cc225023 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -116,7 +116,7 @@ public:
   //! \param[in] r Position to get indices for
   //! \param[out] ijk Array of mesh indices
   //! \param[out] in_mesh Whether position is in mesh
-  virtual void get_indices(Position r, int* ijk, bool* in_mesh) const = 0;
+  virtual void get_indices(Position r, int* ijk, bool* in_mesh) const;
 
   //! Get mesh indices corresponding to a mesh bin
   //
@@ -124,12 +124,19 @@ public:
   //! \param[out] ijk Mesh indices
   virtual void get_indices_from_bin(int bin, int* ijk) const = 0;
 
+  //! Get mesh index in a particular direction
+  //!
+  //! \param[in] r Position to get index for
+  //! \param[in] i Direction index
+  virtual int get_index_in_direction(Position r, int i) const = 0;
+
   //! Get a label for the mesh bin
   std::string bin_label(int bin) const override;
 
   // Data members
   xt::xtensor<double, 1> lower_left_; //!< Lower-left coordinates of mesh
   xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
+  xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
 };
 
 //==============================================================================
@@ -155,10 +162,10 @@ public:
 
   int get_bin_from_indices(const int* ijk) const override;
 
-  void get_indices(Position r, int* ijk, bool* in_mesh) const override;
-
   void get_indices_from_bin(int bin, int* ijk) const override;
 
+  int get_index_in_direction(Position r, int i) const override;
+
   int n_bins() const override;
 
   int n_surface_bins() const override;
@@ -189,7 +196,6 @@ public:
   // Data members
 
   double volume_frac_; //!< Volume fraction of each mesh element
-  xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
   xt::xtensor<double, 1> width_; //!< Width of each mesh element
 
 private:
@@ -217,10 +223,10 @@ public:
 
   int get_bin_from_indices(const int* ijk) const override;
 
-  void get_indices(Position r, int* ijk, bool* in_mesh) const override;
-
   void get_indices_from_bin(int bin, int* ijk) const override;
 
+  int get_index_in_direction(Position r, int i) const override;
+
   int n_bins() const override;
 
   int n_surface_bins() const override;
@@ -240,9 +246,6 @@ public:
   //! \return Whether the line segment connecting r0 and r1 intersects mesh
   bool intersects(Position& r0, Position r1, int* ijk) const;
 
-  // Data members
-  xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
-
 private:
   std::vector<std::vector<double>> grid_;
 };
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 5122c2b617..5c4da64617 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -102,6 +102,17 @@ StructuredMesh::bin_label(int bin) const {
   }
 }
 
+void
+StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
+{
+  *in_mesh = true;
+  for (int i = 0; i < n_dimension_; ++i) {
+    ijk[i] = get_index_in_direction(r, i);
+
+    if (ijk[i] < 1 || ijk[i] > shape_[i]) *in_mesh = false;
+  }
+}
+
 //==============================================================================
 // RegularMesh implementation
 //==============================================================================
@@ -213,16 +224,9 @@ int RegularMesh::get_bin_from_indices(const int* ijk) const
   }
 }
 
-void RegularMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
+int RegularMesh::get_index_in_direction(Position r, int i) const
 {
-  // Find particle in mesh
-  *in_mesh = true;
-  for (int i = 0; i < n_dimension_; ++i) {
-    ijk[i] = std::ceil((r[i] - lower_left_[i]) / width_[i]);
-
-    // Check if indices are within bounds
-    if (ijk[i] < 1 || ijk[i] > shape_[i]) *in_mesh = false;
-  }
+  return std::ceil((r[i] - lower_left_[i]) / width_[i]);
 }
 
 void RegularMesh::get_indices_from_bin(int bin, int* ijk) const
@@ -1175,18 +1179,9 @@ int RectilinearMesh::get_bin_from_indices(const int* ijk) const
   return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0] - 1;
 }
 
-void RectilinearMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
+int RectilinearMesh::get_index_in_direction(Position r, int i) const
 {
-  *in_mesh = true;
-
-  for (int i = 0; i < 3; ++i) {
-    if (r[i] < grid_[i].front() || r[i] > grid_[i].back()) {
-      ijk[i] = -1;
-      *in_mesh = false;
-    } else {
-      ijk[i] = lower_bound_index(grid_[i].begin(), grid_[i].end(), r[i]) + 1;
-    }
-  }
+  return lower_bound_index(grid_[i].begin(), grid_[i].end(), r[i]) + 1;
 }
 
 void RectilinearMesh::get_indices_from_bin(int bin, int* ijk) const

From 7d0838d27b11def6adc2c627387602e777b234c0 Mon Sep 17 00:00:00 2001
From: aprilnovak <novak@berkeley.edu>
Date: Thu, 15 Oct 2020 16:22:27 -0500
Subject: [PATCH 2/8] Move get_bin_from_indices into StructuredMesh base class.
 Refs #1695

---
 include/openmc/mesh.h |  6 +-----
 src/mesh.cpp          | 36 +++++++++++++++---------------------
 2 files changed, 16 insertions(+), 26 deletions(-)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 75cc225023..cfcc137e8e 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -109,7 +109,7 @@ public:
   //
   //! \param[in] Array of mesh indices
   //! \return Mesh bin
-  virtual int get_bin_from_indices(const int* ijk) const = 0;
+  virtual int get_bin_from_indices(const int* ijk) const;
 
   //! Get mesh indices given a position
   //
@@ -160,8 +160,6 @@ public:
 
   int get_bin(Position r) const override;
 
-  int get_bin_from_indices(const int* ijk) const override;
-
   void get_indices_from_bin(int bin, int* ijk) const override;
 
   int get_index_in_direction(Position r, int i) const override;
@@ -221,8 +219,6 @@ public:
 
   int get_bin(Position r) const override;
 
-  int get_bin_from_indices(const int* ijk) const override;
-
   void get_indices_from_bin(int bin, int* ijk) const override;
 
   int get_index_in_direction(Position r, int i) const override;
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 5c4da64617..fd5b30466d 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -102,8 +102,7 @@ StructuredMesh::bin_label(int bin) const {
   }
 }
 
-void
-StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
+void StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
 {
   *in_mesh = true;
   for (int i = 0; i < n_dimension_; ++i) {
@@ -113,6 +112,20 @@ StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
   }
 }
 
+int StructuredMesh::get_bin_from_indices(const int* ijk) const
+{
+  switch (n_dimension_) {
+  case 1:
+    return ijk[0] - 1;
+  case 2:
+    return (ijk[1] - 1)*shape_[0] + ijk[0] - 1;
+  case 3:
+    return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0] - 1;
+  default:
+    throw std::runtime_error{"Invalid number of mesh dimensions"};
+  }
+}
+
 //==============================================================================
 // RegularMesh implementation
 //==============================================================================
@@ -210,20 +223,6 @@ int RegularMesh::get_bin(Position r) const
   return get_bin_from_indices(ijk.data());
 }
 
-int RegularMesh::get_bin_from_indices(const int* ijk) const
-{
-  switch (n_dimension_) {
-  case 1:
-    return ijk[0] - 1;
-  case 2:
-    return (ijk[1] - 1)*shape_[0] + ijk[0] - 1;
-  case 3:
-    return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0] - 1;
-  default:
-    throw std::runtime_error{"Invalid number of mesh dimensions"};
-  }
-}
-
 int RegularMesh::get_index_in_direction(Position r, int i) const
 {
   return std::ceil((r[i] - lower_left_[i]) / width_[i]);
@@ -1174,11 +1173,6 @@ int RectilinearMesh::get_bin(Position r) const
   return get_bin_from_indices(ijk);
 }
 
-int RectilinearMesh::get_bin_from_indices(const int* ijk) const
-{
-  return ((ijk[2] - 1)*shape_[1] + (ijk[1] - 1))*shape_[0] + ijk[0] - 1;
-}
-
 int RectilinearMesh::get_index_in_direction(Position r, int i) const
 {
   return lower_bound_index(grid_[i].begin(), grid_[i].end(), r[i]) + 1;

From 45b7dfb5a03b0e48620a7c1ac376e2dce57fca84 Mon Sep 17 00:00:00 2001
From: aprilnovak <novak@berkeley.edu>
Date: Thu, 15 Oct 2020 16:26:36 -0500
Subject: [PATCH 3/8] Move implementation of get_indices_from_bin to
 StructuredMesh base class. Refs #1695

---
 include/openmc/mesh.h |  6 +-----
 src/mesh.cpp          | 35 ++++++++++++++---------------------
 2 files changed, 15 insertions(+), 26 deletions(-)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index cfcc137e8e..9f1b80c035 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -122,7 +122,7 @@ public:
   //
   //! \param[in] bin Mesh bin
   //! \param[out] ijk Mesh indices
-  virtual void get_indices_from_bin(int bin, int* ijk) const = 0;
+  virtual void get_indices_from_bin(int bin, int* ijk) const;
 
   //! Get mesh index in a particular direction
   //!
@@ -160,8 +160,6 @@ public:
 
   int get_bin(Position r) const override;
 
-  void get_indices_from_bin(int bin, int* ijk) const override;
-
   int get_index_in_direction(Position r, int i) const override;
 
   int n_bins() const override;
@@ -219,8 +217,6 @@ public:
 
   int get_bin(Position r) const override;
 
-  void get_indices_from_bin(int bin, int* ijk) const override;
-
   int get_index_in_direction(Position r, int i) const override;
 
   int n_bins() const override;
diff --git a/src/mesh.cpp b/src/mesh.cpp
index fd5b30466d..17a87c47df 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -126,6 +126,20 @@ int StructuredMesh::get_bin_from_indices(const int* ijk) const
   }
 }
 
+void StructuredMesh::get_indices_from_bin(int bin, int* ijk) const
+{
+  if (n_dimension_ == 1) {
+    ijk[0] = bin + 1;
+  } else if (n_dimension_ == 2) {
+    ijk[0] = bin % shape_[0] + 1;
+    ijk[1] = bin / shape_[0] + 1;
+  } else if (n_dimension_ == 3) {
+    ijk[0] = bin % shape_[0] + 1;
+    ijk[1] = (bin % (shape_[0] * shape_[1])) / shape_[0] + 1;
+    ijk[2] = bin / (shape_[0] * shape_[1]) + 1;
+  }
+}
+
 //==============================================================================
 // RegularMesh implementation
 //==============================================================================
@@ -228,20 +242,6 @@ int RegularMesh::get_index_in_direction(Position r, int i) const
   return std::ceil((r[i] - lower_left_[i]) / width_[i]);
 }
 
-void RegularMesh::get_indices_from_bin(int bin, int* ijk) const
-{
-  if (n_dimension_ == 1) {
-    ijk[0] = bin + 1;
-  } else if (n_dimension_ == 2) {
-    ijk[0] = bin % shape_[0] + 1;
-    ijk[1] = bin / shape_[0] + 1;
-  } else if (n_dimension_ == 3) {
-    ijk[0] = bin % shape_[0] + 1;
-    ijk[1] = (bin % (shape_[0] * shape_[1])) / shape_[0] + 1;
-    ijk[2] = bin / (shape_[0] * shape_[1]) + 1;
-  }
-}
-
 int RegularMesh::n_bins() const
 {
   int n_bins = 1;
@@ -1178,13 +1178,6 @@ int RectilinearMesh::get_index_in_direction(Position r, int i) const
   return lower_bound_index(grid_[i].begin(), grid_[i].end(), r[i]) + 1;
 }
 
-void RectilinearMesh::get_indices_from_bin(int bin, int* ijk) const
-{
-  ijk[0] = bin % shape_[0] + 1;
-  ijk[1] = (bin % (shape_[0] * shape_[1])) / shape_[0] + 1;
-  ijk[2] = bin / (shape_[0] * shape_[1]) + 1;
-}
-
 int RectilinearMesh::n_bins() const
 {
   return xt::prod(shape_)();

From bcc072e0a8774f2cda8e3708ec73cba6631615d4 Mon Sep 17 00:00:00 2001
From: aprilnovak <novak@berkeley.edu>
Date: Thu, 15 Oct 2020 16:57:23 -0500
Subject: [PATCH 4/8] Move intersects method up to base class and combine the
 functionality using new get_index_in_direction function. Refs #1695

---
 include/openmc/mesh.h |  40 ++-
 src/mesh.cpp          | 608 +++++++++++++++++-------------------------
 2 files changed, 264 insertions(+), 384 deletions(-)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 9f1b80c035..e619c43705 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -126,9 +126,17 @@ public:
 
   //! Get mesh index in a particular direction
   //!
-  //! \param[in] r Position to get index for
+  //! \param[in] r Coordinate to get index for
   //! \param[in] i Direction index
-  virtual int get_index_in_direction(Position r, int i) const = 0;
+  virtual int get_index_in_direction(double r, int i) const = 0;
+
+  //! Check where a line segment intersects the mesh and if it intersects at all
+  //
+  //! \param[in,out] r0 In: starting position, out: intersection point
+  //! \param[in] r1 Ending position
+  //! \param[out] ijk Indices of the mesh bin containing the intersection point
+  //! \return Whether the line segment connecting r0 and r1 intersects mesh
+  virtual bool intersects(Position& r0, Position r1, int* ijk) const;
 
   //! Get a label for the mesh bin
   std::string bin_label(int bin) const override;
@@ -137,6 +145,11 @@ public:
   xt::xtensor<double, 1> lower_left_; //!< Lower-left coordinates of mesh
   xt::xtensor<double, 1> upper_right_; //!< Upper-right coordinates of mesh
   xt::xtensor<int, 1> shape_; //!< Number of mesh elements in each dimension
+
+protected:
+  virtual bool intersects_1d(Position& r0, Position r1, int* ijk) const;
+  virtual bool intersects_2d(Position& r0, Position r1, int* ijk) const;
+  virtual bool intersects_3d(Position& r0, Position r1, int* ijk) const;
 };
 
 //==============================================================================
@@ -173,14 +186,6 @@ public:
 
   // New methods
 
-  //! Check where a line segment intersects the mesh and if it intersects at all
-  //
-  //! \param[in,out] r0 In: starting position, out: intersection point
-  //! \param[in] r1 Ending position
-  //! \param[out] ijk Indices of the mesh bin containing the intersection point
-  //! \return Whether the line segment connecting r0 and r1 intersects mesh
-  bool intersects(Position& r0, Position r1, int* ijk) const;
-
   //! Count number of bank sites in each mesh bin / energy bin
   //
   //! \param[in] bank Array of bank sites
@@ -193,11 +198,6 @@ public:
 
   double volume_frac_; //!< Volume fraction of each mesh element
   xt::xtensor<double, 1> width_; //!< Width of each mesh element
-
-private:
-  bool intersects_1d(Position& r0, Position r1, int* ijk) const;
-  bool intersects_2d(Position& r0, Position r1, int* ijk) const;
-  bool intersects_3d(Position& r0, Position r1, int* ijk) const;
 };
 
 
@@ -228,16 +228,6 @@ public:
 
   void to_hdf5(hid_t group) const override;
 
-  // New methods
-
-  //! Check where a line segment intersects the mesh and if it intersects at all
-  //
-  //! \param[in,out] r0 In: starting position, out: intersection point
-  //! \param[in] r1 Ending position
-  //! \param[out] ijk Indices of the mesh bin containing the intersection point
-  //! \return Whether the line segment connecting r0 and r1 intersects mesh
-  bool intersects(Position& r0, Position r1, int* ijk) const;
-
 private:
   std::vector<std::vector<double>> grid_;
 };
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 17a87c47df..c4158df233 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -106,7 +106,7 @@ void StructuredMesh::get_indices(Position r, int* ijk, bool* in_mesh) const
 {
   *in_mesh = true;
   for (int i = 0; i < n_dimension_; ++i) {
-    ijk[i] = get_index_in_direction(r, i);
+    ijk[i] = get_index_in_direction(r[i], i);
 
     if (ijk[i] < 1 || ijk[i] > shape_[i]) *in_mesh = false;
   }
@@ -140,6 +140,250 @@ void StructuredMesh::get_indices_from_bin(int bin, int* ijk) const
   }
 }
 
+bool StructuredMesh::intersects(Position& r0, Position r1, int* ijk) const
+{
+  switch(n_dimension_) {
+  case 1:
+    return intersects_1d(r0, r1, ijk);
+  case 2:
+    return intersects_2d(r0, r1, ijk);
+  case 3:
+    return intersects_3d(r0, r1, ijk);
+  default:
+    throw std::runtime_error{"Invalid number of mesh dimensions."};
+  }
+}
+
+bool StructuredMesh::intersects_1d(Position& r0, Position r1, int* ijk) const
+{
+  // Copy coordinates of starting point
+  double x0 = r0.x;
+  double y0 = r0.y;
+  double z0 = r0.z;
+
+  // Copy coordinates of ending point
+  double x1 = r1.x;
+  double y1 = r1.y;
+  double z1 = r1.z;
+
+  // Copy coordinates of mesh lower_left and upper_right
+  double xm0 = lower_left_[0];
+  double xm1 = upper_right_[0];
+
+  double min_dist = INFTY;
+
+  // Check if line intersects left surface -- calculate the intersection point
+  // (y,z)
+  if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
+    double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
+    double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
+    if (check_intersection_point(xm0, x0, yi, yi, zi, zi, r0, min_dist)) {
+      ijk[0] = 1;
+    }
+  }
+
+  // Check if line intersects right surface -- calculate the intersection point
+  // (y,z)
+  if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
+    double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
+    double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
+    if (check_intersection_point(xm1, x0, yi, yi, zi, zi, r0, min_dist)) {
+      ijk[0] = shape_[0];
+    }
+  }
+
+  return min_dist < INFTY;
+}
+
+bool StructuredMesh::intersects_2d(Position& r0, Position r1, int* ijk) const
+{
+  // Copy coordinates of starting point
+  double x0 = r0.x;
+  double y0 = r0.y;
+  double z0 = r0.z;
+
+  // Copy coordinates of ending point
+  double x1 = r1.x;
+  double y1 = r1.y;
+  double z1 = r1.z;
+
+  // Copy coordinates of mesh lower_left
+  double xm0 = lower_left_[0];
+  double ym0 = lower_left_[1];
+
+  // Copy coordinates of mesh upper_right
+  double xm1 = upper_right_[0];
+  double ym1 = upper_right_[1];
+
+  double min_dist = INFTY;
+
+  // Check if line intersects left surface -- calculate the intersection point
+  // (y,z)
+  if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
+    double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
+    double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
+    if (yi >= ym0 && yi < ym1) {
+      if (check_intersection_point(xm0, x0, yi, y0, zi, zi, r0, min_dist)) {
+        ijk[0] = 1;
+        ijk[1] = get_index_in_direction(yi, 1);
+      }
+    }
+  }
+
+  // Check if line intersects back surface -- calculate the intersection point
+  // (x,z)
+  if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
+    double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
+    double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
+    if (xi >= xm0 && xi < xm1) {
+      if (check_intersection_point(xi, x0, ym0, y0, zi, zi, r0, min_dist)) {
+        ijk[0] = get_index_in_direction(xi, 0);
+        ijk[1] = 1;
+      }
+    }
+  }
+
+  // Check if line intersects right surface -- calculate the intersection point
+  // (y,z)
+  if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
+    double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
+    double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
+    if (yi >= ym0 && yi < ym1) {
+      if (check_intersection_point(xm1, x0, yi, y0, zi, zi, r0, min_dist)) {
+        ijk[0] = shape_[0];
+        ijk[1] = get_index_in_direction(yi, 1);
+      }
+    }
+  }
+
+  // Check if line intersects front surface -- calculate the intersection point
+  // (x,z)
+  if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
+    double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
+    double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
+    if (xi >= xm0 && xi < xm1) {
+      if (check_intersection_point(xi, x0, ym1, y0, zi, zi, r0, min_dist)) {
+        ijk[0] = get_index_in_direction(xi, 0);
+        ijk[1] = shape_[1];
+      }
+    }
+  }
+
+  return min_dist < INFTY;
+}
+
+bool StructuredMesh::intersects_3d(Position& r0, Position r1, int* ijk) const
+{
+  // Copy coordinates of starting point
+  double x0 = r0.x;
+  double y0 = r0.y;
+  double z0 = r0.z;
+
+  // Copy coordinates of ending point
+  double x1 = r1.x;
+  double y1 = r1.y;
+  double z1 = r1.z;
+
+  // Copy coordinates of mesh lower_left
+  double xm0 = lower_left_[0];
+  double ym0 = lower_left_[1];
+  double zm0 = lower_left_[2];
+
+  // Copy coordinates of mesh upper_right
+  double xm1 = upper_right_[0];
+  double ym1 = upper_right_[1];
+  double zm1 = upper_right_[2];
+
+  double min_dist = INFTY;
+
+  // Check if line intersects left surface -- calculate the intersection point
+  // (y,z)
+  if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
+    double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
+    double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
+    if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
+      if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) {
+        ijk[0] = 1;
+        ijk[1] = get_index_in_direction(yi, 1);
+        ijk[2] = get_index_in_direction(zi, 2);
+      }
+    }
+  }
+
+  // Check if line intersects back surface -- calculate the intersection point
+  // (x,z)
+  if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
+    double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
+    double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
+    if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
+      if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) {
+        ijk[0] = get_index_in_direction(xi, 0);
+        ijk[1] = 1;
+        ijk[2] = get_index_in_direction(zi, 2);
+      }
+    }
+  }
+
+  // Check if line intersects bottom surface -- calculate the intersection
+  // point (x,y)
+  if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) {
+    double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0);
+    double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0);
+    if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
+      if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) {
+        ijk[0] = get_index_in_direction(xi, 0);
+        ijk[1] = get_index_in_direction(yi, 1);
+        ijk[2] = 1;
+      }
+    }
+  }
+
+  // Check if line intersects right surface -- calculate the intersection point
+  // (y,z)
+  if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
+    double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
+    double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
+    if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
+      if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) {
+        ijk[0] = shape_[0];
+        ijk[1] = get_index_in_direction(yi, 1);
+        ijk[2] = get_index_in_direction(zi, 2);
+      }
+    }
+  }
+
+  // Check if line intersects front surface -- calculate the intersection point
+  // (x,z)
+  if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
+    double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
+    double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
+    if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
+      if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) {
+        ijk[0] = get_index_in_direction(xi, 0);
+        ijk[1] = shape_[1];
+        ijk[2] = get_index_in_direction(zi, 2);
+      }
+    }
+  }
+
+  // Check if line intersects top surface -- calculate the intersection point
+  // (x,y)
+  if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) {
+    double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0);
+    double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0);
+    if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
+      if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) {
+        ijk[0] = get_index_in_direction(xi, 0);
+        ijk[1] = get_index_in_direction(yi, 1);
+        ijk[2] = shape_[2];
+      }
+    }
+  }
+
+  return min_dist < INFTY;
+}
+
+
 //==============================================================================
 // RegularMesh implementation
 //==============================================================================
@@ -237,9 +481,9 @@ int RegularMesh::get_bin(Position r) const
   return get_bin_from_indices(ijk.data());
 }
 
-int RegularMesh::get_index_in_direction(Position r, int i) const
+int RegularMesh::get_index_in_direction(double r, int i) const
 {
-  return std::ceil((r[i] - lower_left_[i]) / width_[i]);
+  return std::ceil((r - lower_left_[i]) / width_[i]);
 }
 
 int RegularMesh::n_bins() const
@@ -254,249 +498,6 @@ int RegularMesh::n_surface_bins() const
   return 4 * n_dimension_ * n_bins();
 }
 
-bool RegularMesh::intersects(Position& r0, Position r1, int* ijk) const
-{
-  switch(n_dimension_) {
-  case 1:
-    return intersects_1d(r0, r1, ijk);
-  case 2:
-    return intersects_2d(r0, r1, ijk);
-  case 3:
-    return intersects_3d(r0, r1, ijk);
-  default:
-    throw std::runtime_error{"Invalid number of mesh dimensions."};
-  }
-}
-
-bool RegularMesh::intersects_1d(Position& r0, Position r1, int* ijk) const
-{
-  // Copy coordinates of starting point
-  double x0 = r0.x;
-  double y0 = r0.y;
-  double z0 = r0.z;
-
-  // Copy coordinates of ending point
-  double x1 = r1.x;
-  double y1 = r1.y;
-  double z1 = r1.z;
-
-  // Copy coordinates of mesh lower_left and upper_right
-  double xm0 = lower_left_[0];
-  double xm1 = upper_right_[0];
-
-  double min_dist = INFTY;
-
-  // Check if line intersects left surface -- calculate the intersection point
-  // (y,z)
-  if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
-    double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
-    double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
-    if (check_intersection_point(xm0, x0, yi, yi, zi, zi, r0, min_dist)) {
-      ijk[0] = 1;
-    }
-  }
-
-  // Check if line intersects right surface -- calculate the intersection point
-  // (y,z)
-  if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
-    double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
-    double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
-    if (check_intersection_point(xm1, x0, yi, yi, zi, zi, r0, min_dist)) {
-      ijk[0] = shape_[0];
-    }
-  }
-
-  return min_dist < INFTY;
-}
-
-bool RegularMesh::intersects_2d(Position& r0, Position r1, int* ijk) const
-{
-  // Copy coordinates of starting point
-  double x0 = r0.x;
-  double y0 = r0.y;
-  double z0 = r0.z;
-
-  // Copy coordinates of ending point
-  double x1 = r1.x;
-  double y1 = r1.y;
-  double z1 = r1.z;
-
-  // Copy coordinates of mesh lower_left
-  double xm0 = lower_left_[0];
-  double ym0 = lower_left_[1];
-
-  // Copy coordinates of mesh upper_right
-  double xm1 = upper_right_[0];
-  double ym1 = upper_right_[1];
-
-  double min_dist = INFTY;
-
-  // Check if line intersects left surface -- calculate the intersection point
-  // (y,z)
-  if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
-    double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
-    double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
-    if (yi >= ym0 && yi < ym1) {
-      if (check_intersection_point(xm0, x0, yi, y0, zi, zi, r0, min_dist)) {
-        ijk[0] = 1;
-        ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
-      }
-    }
-  }
-
-  // Check if line intersects back surface -- calculate the intersection point
-  // (x,z)
-  if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
-    double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
-    double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
-    if (xi >= xm0 && xi < xm1) {
-      if (check_intersection_point(xi, x0, ym0, y0, zi, zi, r0, min_dist)) {
-        ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
-        ijk[1] = 1;
-      }
-    }
-  }
-
-  // Check if line intersects right surface -- calculate the intersection point
-  // (y,z)
-  if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
-    double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
-    double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
-    if (yi >= ym0 && yi < ym1) {
-      if (check_intersection_point(xm1, x0, yi, y0, zi, zi, r0, min_dist)) {
-        ijk[0] = shape_[0];
-        ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
-      }
-    }
-  }
-
-  // Check if line intersects front surface -- calculate the intersection point
-  // (x,z)
-  if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
-    double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
-    double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
-    if (xi >= xm0 && xi < xm1) {
-      if (check_intersection_point(xi, x0, ym1, y0, zi, zi, r0, min_dist)) {
-        ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
-        ijk[1] = shape_[1];
-      }
-    }
-  }
-
-  return min_dist < INFTY;
-}
-
-bool RegularMesh::intersects_3d(Position& r0, Position r1, int* ijk) const
-{
-  // Copy coordinates of starting point
-  double x0 = r0.x;
-  double y0 = r0.y;
-  double z0 = r0.z;
-
-  // Copy coordinates of ending point
-  double x1 = r1.x;
-  double y1 = r1.y;
-  double z1 = r1.z;
-
-  // Copy coordinates of mesh lower_left
-  double xm0 = lower_left_[0];
-  double ym0 = lower_left_[1];
-  double zm0 = lower_left_[2];
-
-  // Copy coordinates of mesh upper_right
-  double xm1 = upper_right_[0];
-  double ym1 = upper_right_[1];
-  double zm1 = upper_right_[2];
-
-  double min_dist = INFTY;
-
-  // Check if line intersects left surface -- calculate the intersection point
-  // (y,z)
-  if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
-    double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
-    double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
-    if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
-      if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) {
-        ijk[0] = 1;
-        ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
-        ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
-      }
-    }
-  }
-
-  // Check if line intersects back surface -- calculate the intersection point
-  // (x,z)
-  if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
-    double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
-    double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
-    if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
-      if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) {
-        ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
-        ijk[1] = 1;
-        ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
-      }
-    }
-  }
-
-  // Check if line intersects bottom surface -- calculate the intersection
-  // point (x,y)
-  if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) {
-    double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0);
-    double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0);
-    if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
-      if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) {
-        ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
-        ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
-        ijk[2] = 1;
-      }
-    }
-  }
-
-  // Check if line intersects right surface -- calculate the intersection point
-  // (y,z)
-  if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
-    double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
-    double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
-    if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
-      if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) {
-        ijk[0] = shape_[0];
-        ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
-        ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
-      }
-    }
-  }
-
-  // Check if line intersects front surface -- calculate the intersection point
-  // (x,z)
-  if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
-    double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
-    double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
-    if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
-      if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) {
-        ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
-        ijk[1] = shape_[1];
-        ijk[2] = std::ceil((zi - lower_left_[2]) / width_[2]);
-      }
-    }
-  }
-
-  // Check if line intersects top surface -- calculate the intersection point
-  // (x,y)
-  if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) {
-    double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0);
-    double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0);
-    if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
-      if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) {
-        ijk[0] = std::ceil((xi - lower_left_[0]) / width_[0]);
-        ijk[1] = std::ceil((yi - lower_left_[1]) / width_[1]);
-        ijk[2] = shape_[2];
-      }
-    }
-  }
-
-  return min_dist < INFTY;
-}
-
 void RegularMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
                                std::vector<double>& lengths) const
 {
@@ -1173,9 +1174,9 @@ int RectilinearMesh::get_bin(Position r) const
   return get_bin_from_indices(ijk);
 }
 
-int RectilinearMesh::get_index_in_direction(Position r, int i) const
+int RectilinearMesh::get_index_in_direction(double r, int i) const
 {
-  return lower_bound_index(grid_[i].begin(), grid_[i].end(), r[i]) + 1;
+  return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
 }
 
 int RectilinearMesh::n_bins() const
@@ -1230,117 +1231,6 @@ void RectilinearMesh::to_hdf5(hid_t group) const
   close_group(mesh_group);
 }
 
-bool RectilinearMesh::intersects(Position& r0, Position r1, int* ijk) const
-{
-  // Copy coordinates of starting point
-  double x0 = r0.x;
-  double y0 = r0.y;
-  double z0 = r0.z;
-
-  // Copy coordinates of ending point
-  double x1 = r1.x;
-  double y1 = r1.y;
-  double z1 = r1.z;
-
-  // Copy coordinates of mesh lower_left
-  double xm0 = grid_[0].front();
-  double ym0 = grid_[1].front();
-  double zm0 = grid_[2].front();
-
-  // Copy coordinates of mesh upper_right
-  double xm1 = grid_[0].back();
-  double ym1 = grid_[1].back();
-  double zm1 = grid_[2].back();
-
-  double min_dist = INFTY;
-
-  // Check if line intersects left surface -- calculate the intersection point
-  // (y,z)
-  if ((x0 < xm0 && x1 > xm0) || (x0 > xm0 && x1 < xm0)) {
-    double yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0);
-    double zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0);
-    if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
-      if (check_intersection_point(xm0, x0, yi, y0, zi, z0, r0, min_dist)) {
-        ijk[0] = 1;
-        ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
-        ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
-      }
-    }
-  }
-
-  // Check if line intersects back surface -- calculate the intersection point
-  // (x,z)
-  if ((y0 < ym0 && y1 > ym0) || (y0 > ym0 && y1 < ym0)) {
-    double xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0);
-    double zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0);
-    if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
-      if (check_intersection_point(xi, x0, ym0, y0, zi, z0, r0, min_dist)) {
-        ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
-        ijk[1] = 1;
-        ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
-      }
-    }
-  }
-
-  // Check if line intersects bottom surface -- calculate the intersection
-  // point (x,y)
-  if ((z0 < zm0 && z1 > zm0) || (z0 > zm0 && z1 < zm0)) {
-    double xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0);
-    double yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0);
-    if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
-      if (check_intersection_point(xi, x0, yi, y0, zm0, z0, r0, min_dist)) {
-        ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
-        ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
-        ijk[2] = 1;
-      }
-    }
-  }
-
-  // Check if line intersects right surface -- calculate the intersection point
-  // (y,z)
-  if ((x0 < xm1 && x1 > xm1) || (x0 > xm1 && x1 < xm1)) {
-    double yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0);
-    double zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0);
-    if (yi >= ym0 && yi < ym1 && zi >= zm0 && zi < zm1) {
-      if (check_intersection_point(xm1, x0, yi, y0, zi, z0, r0, min_dist)) {
-        ijk[0] = shape_[0];
-        ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
-        ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
-      }
-    }
-  }
-
-  // Check if line intersects front surface -- calculate the intersection point
-  // (x,z)
-  if ((y0 < ym1 && y1 > ym1) || (y0 > ym1 && y1 < ym1)) {
-    double xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0);
-    double zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0);
-    if (xi >= xm0 && xi < xm1 && zi >= zm0 && zi < zm1) {
-      if (check_intersection_point(xi, x0, ym1, y0, zi, z0, r0, min_dist)) {
-        ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
-        ijk[1] = shape_[1];
-        ijk[2] = lower_bound_index(grid_[2].begin(), grid_[2].end(), zi) + 1;
-      }
-    }
-  }
-
-  // Check if line intersects top surface -- calculate the intersection point
-  // (x,y)
-  if ((z0 < zm1 && z1 > zm1) || (z0 > zm1 && z1 < zm1)) {
-    double xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0);
-    double yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0);
-    if (xi >= xm0 && xi < xm1 && yi >= ym0 && yi < ym1) {
-      if (check_intersection_point(xi, x0, yi, y0, zm1, z0, r0, min_dist)) {
-        ijk[0] = lower_bound_index(grid_[0].begin(), grid_[0].end(), xi) + 1;
-        ijk[1] = lower_bound_index(grid_[1].begin(), grid_[1].end(), yi) + 1;
-        ijk[2] = shape_[2];
-      }
-    }
-  }
-
-  return min_dist < INFTY;
-}
-
 //==============================================================================
 // Helper functions for the C API
 //==============================================================================

From 9a36ea405aeeed940152882199219cdf6a35207c Mon Sep 17 00:00:00 2001
From: aprilnovak <novak@berkeley.edu>
Date: Thu, 15 Oct 2020 17:14:32 -0500
Subject: [PATCH 5/8] Move get_bin implementation up to StructuredMesh base
 class. Refs #1695

---
 include/openmc/mesh.h | 10 ++++------
 src/mesh.cpp          | 36 ++++++++++++------------------------
 2 files changed, 16 insertions(+), 30 deletions(-)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index e619c43705..98f0e7a83c 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -105,6 +105,8 @@ public:
   StructuredMesh(pugi::xml_node node) : Mesh {node} {};
   virtual ~StructuredMesh() = default;
 
+  int get_bin(Position r) const override;
+
   //! Get bin given mesh indices
   //
   //! \param[in] Array of mesh indices
@@ -171,9 +173,7 @@ public:
   void surface_bins_crossed(const Particle& p, std::vector<int>& bins)
   const override;
 
-  int get_bin(Position r) const override;
-
-  int get_index_in_direction(Position r, int i) const override;
+  int get_index_in_direction(double r, int i) const override;
 
   int n_bins() const override;
 
@@ -215,9 +215,7 @@ public:
   void surface_bins_crossed(const Particle& p, std::vector<int>& bins)
   const override;
 
-  int get_bin(Position r) const override;
-
-  int get_index_in_direction(Position r, int i) const override;
+  int get_index_in_direction(double r, int i) const override;
 
   int n_bins() const override;
 
diff --git a/src/mesh.cpp b/src/mesh.cpp
index c4158df233..68a116e3af 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -140,6 +140,18 @@ void StructuredMesh::get_indices_from_bin(int bin, int* ijk) const
   }
 }
 
+int StructuredMesh::get_bin(Position r) const
+{
+  // Determine indices
+  std::vector<int> ijk(n_dimension_);
+  bool in_mesh;
+  get_indices(r, ijk.data(), &in_mesh);
+  if (!in_mesh) return -1;
+
+  // Convert indices to bin
+  return get_bin_from_indices(ijk.data());
+}
+
 bool StructuredMesh::intersects(Position& r0, Position r1, int* ijk) const
 {
   switch(n_dimension_) {
@@ -469,18 +481,6 @@ RegularMesh::RegularMesh(pugi::xml_node node)
   volume_frac_ = 1.0/xt::prod(shape_)();
 }
 
-int RegularMesh::get_bin(Position r) const
-{
-  // Determine indices
-  std::vector<int> ijk(n_dimension_);
-  bool in_mesh;
-  get_indices(r, ijk.data(), &in_mesh);
-  if (!in_mesh) return -1;
-
-  // Convert indices to bin
-  return get_bin_from_indices(ijk.data());
-}
-
 int RegularMesh::get_index_in_direction(double r, int i) const
 {
   return std::ceil((r - lower_left_[i]) / width_[i]);
@@ -1162,18 +1162,6 @@ void RectilinearMesh::surface_bins_crossed(const Particle& p,
   }
 }
 
-int RectilinearMesh::get_bin(Position r) const
-{
-  // Determine indices
-  int ijk[3];
-  bool in_mesh;
-  get_indices(r, ijk, &in_mesh);
-  if (!in_mesh) return -1;
-
-  // Convert indices to bin
-  return get_bin_from_indices(ijk);
-}
-
 int RectilinearMesh::get_index_in_direction(double r, int i) const
 {
   return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;

From 008a73fa361308987e7141dcdb98be327bbc5b2d Mon Sep 17 00:00:00 2001
From: aprilnovak <novak@berkeley.edu>
Date: Thu, 15 Oct 2020 17:18:32 -0500
Subject: [PATCH 6/8] Move n_bins and n_surface_bins implementations to
 StructuredMesh base class. Refs #1695

---
 include/openmc/mesh.h | 12 ++++--------
 src/mesh.cpp          | 32 ++++++++++----------------------
 2 files changed, 14 insertions(+), 30 deletions(-)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 98f0e7a83c..f3f84ac6b7 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -107,6 +107,10 @@ public:
 
   int get_bin(Position r) const override;
 
+  int n_bins() const override;
+
+  int n_surface_bins() const override;
+
   //! Get bin given mesh indices
   //
   //! \param[in] Array of mesh indices
@@ -175,10 +179,6 @@ public:
 
   int get_index_in_direction(double r, int i) const override;
 
-  int n_bins() const override;
-
-  int n_surface_bins() const override;
-
   std::pair<std::vector<double>, std::vector<double>>
   plot(Position plot_ll, Position plot_ur) const override;
 
@@ -217,10 +217,6 @@ public:
 
   int get_index_in_direction(double r, int i) const override;
 
-  int n_bins() const override;
-
-  int n_surface_bins() const override;
-
   std::pair<std::vector<double>, std::vector<double>>
   plot(Position plot_ll, Position plot_ur) const override;
 
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 68a116e3af..ca7eb53ab9 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -152,6 +152,16 @@ int StructuredMesh::get_bin(Position r) const
   return get_bin_from_indices(ijk.data());
 }
 
+int StructuredMesh::n_bins() const
+{
+  return xt::prod(shape_)();
+}
+
+int StructuredMesh::n_surface_bins() const
+{
+  return 4 * n_dimension_ * n_bins();
+}
+
 bool StructuredMesh::intersects(Position& r0, Position r1, int* ijk) const
 {
   switch(n_dimension_) {
@@ -486,18 +496,6 @@ int RegularMesh::get_index_in_direction(double r, int i) const
   return std::ceil((r - lower_left_[i]) / width_[i]);
 }
 
-int RegularMesh::n_bins() const
-{
-  int n_bins = 1;
-  for (auto dim : shape_) n_bins *= dim;
-  return n_bins;
-}
-
-int RegularMesh::n_surface_bins() const
-{
-  return 4 * n_dimension_ * n_bins();
-}
-
 void RegularMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
                                std::vector<double>& lengths) const
 {
@@ -1167,16 +1165,6 @@ int RectilinearMesh::get_index_in_direction(double r, int i) const
   return lower_bound_index(grid_[i].begin(), grid_[i].end(), r) + 1;
 }
 
-int RectilinearMesh::n_bins() const
-{
-  return xt::prod(shape_)();
-}
-
-int RectilinearMesh::n_surface_bins() const
-{
-  return 4 * n_dimension_ * n_bins();
-}
-
 std::pair<std::vector<double>, std::vector<double>>
 RectilinearMesh::plot(Position plot_ll, Position plot_ur) const
 {

From 1207c901be203110b6bc6b0a5b8b791e10a0a1b9 Mon Sep 17 00:00:00 2001
From: aprilnovak <novak@berkeley.edu>
Date: Thu, 15 Oct 2020 18:01:02 -0500
Subject: [PATCH 7/8] Move bins_crossed implementation into StructuredMesh base
 class. Refs #1695

---
 include/openmc/mesh.h |  29 +++-
 src/mesh.cpp          | 318 ++++++++++++++++--------------------------
 2 files changed, 141 insertions(+), 206 deletions(-)

diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index f3f84ac6b7..5a8367b8bb 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -111,6 +111,9 @@ public:
 
   int n_surface_bins() const override;
 
+  void bins_crossed(const Particle& p, std::vector<int>& bins,
+                    std::vector<double>& lengths) const override;
+
   //! Get bin given mesh indices
   //
   //! \param[in] Array of mesh indices
@@ -144,6 +147,18 @@ public:
   //! \return Whether the line segment connecting r0 and r1 intersects mesh
   virtual bool intersects(Position& r0, Position r1, int* ijk) const;
 
+  //! Get the coordinate for the mesh grid boundary in the positive direction
+  //!
+  //! \param[in] ijk Array of mesh indices
+  //! \param[in] i Direction index
+  virtual double positive_grid_boundary(int* ijk, int i) const = 0;
+
+  //! Get the coordinate for the mesh grid boundary in the negative direction
+  //!
+  //! \param[in] ijk Array of mesh indices
+  //! \param[in] i Direction index
+  virtual double negative_grid_boundary(int* ijk, int i) const = 0;
+
   //! Get a label for the mesh bin
   std::string bin_label(int bin) const override;
 
@@ -171,14 +186,15 @@ public:
 
   // Overriden methods
 
-  void bins_crossed(const Particle& p, std::vector<int>& bins,
-                    std::vector<double>& lengths) const override;
-
   void surface_bins_crossed(const Particle& p, std::vector<int>& bins)
   const override;
 
   int get_index_in_direction(double r, int i) const override;
 
+  double positive_grid_boundary(int* ijk, int i) const override;
+
+  double negative_grid_boundary(int* ijk, int i) const override;
+
   std::pair<std::vector<double>, std::vector<double>>
   plot(Position plot_ll, Position plot_ur) const override;
 
@@ -209,14 +225,15 @@ public:
 
   // Overriden methods
 
-  void bins_crossed(const Particle& p, std::vector<int>& bins,
-                    std::vector<double>& lengths) const override;
-
   void surface_bins_crossed(const Particle& p, std::vector<int>& bins)
   const override;
 
   int get_index_in_direction(double r, int i) const override;
 
+  double positive_grid_boundary(int* ijk, int i) const override;
+
+  double negative_grid_boundary(int* ijk, int i) const override;
+
   std::pair<std::vector<double>, std::vector<double>>
   plot(Position plot_ll, Position plot_ur) const override;
 
diff --git a/src/mesh.cpp b/src/mesh.cpp
index ca7eb53ab9..d3b671263a 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -405,99 +405,7 @@ bool StructuredMesh::intersects_3d(Position& r0, Position r1, int* ijk) const
   return min_dist < INFTY;
 }
 
-
-//==============================================================================
-// RegularMesh implementation
-//==============================================================================
-
-RegularMesh::RegularMesh(pugi::xml_node node)
-  : StructuredMesh {node}
-{
-  // Determine number of dimensions for mesh
-  if (check_for_node(node, "dimension")) {
-    shape_ = get_node_xarray<int>(node, "dimension");
-    int n = n_dimension_ = shape_.size();
-    if (n != 1 && n != 2 && n != 3) {
-      fatal_error("Mesh must be one, two, or three dimensions.");
-    }
-
-    // Check that dimensions are all greater than zero
-    if (xt::any(shape_ <= 0)) {
-      fatal_error("All entries on the <dimension> element for a tally "
-        "mesh must be positive.");
-    }
-  } else {
-    fatal_error("Must specify <dimension> on a mesh.");
-  }
-
-  // Check for lower-left coordinates
-  if (check_for_node(node, "lower_left")) {
-    // Read mesh lower-left corner location
-    lower_left_ = get_node_xarray<double>(node, "lower_left");
-
-    // Make sure lower_left and dimension match
-    if (n_dimension_ != lower_left_.size()) {
-      fatal_error("Number of entries on <lower_left> must be the same "
-        "as the number of entries on <dimension>.");
-    }
-  } else {
-    fatal_error("Must specify <lower_left> on a mesh.");
-  }
-
-  if (check_for_node(node, "width")) {
-    // Make sure both upper-right or width were specified
-    if (check_for_node(node, "upper_right")) {
-      fatal_error("Cannot specify both <upper_right> and <width> on a mesh.");
-    }
-
-    width_ = get_node_xarray<double>(node, "width");
-
-    // Check to ensure width has same dimensions
-    if (n_dimension_ != width_.size()) {
-      fatal_error("Number of entries on <width> must be the same as "
-        "the number of entries on <lower_left>.");
-    }
-
-    // Check for negative widths
-    if (xt::any(width_ < 0.0)) {
-      fatal_error("Cannot have a negative <width> on a tally mesh.");
-    }
-
-    // Set width and upper right coordinate
-    upper_right_ = xt::eval(lower_left_ + shape_ * width_);
-
-  } else if (check_for_node(node, "upper_right")) {
-    upper_right_ = get_node_xarray<double>(node, "upper_right");
-
-    // Check to ensure upper right has same dimensions
-    if (n_dimension_ != upper_right_.size()) {
-      fatal_error("Number of entries on <upper_right> must be the "
-        "same as the number of entries on <lower_left>.");
-    }
-
-    // Check that upper-right is above lower-left
-    if (xt::any(upper_right_ < lower_left_)) {
-      fatal_error("The <upper_right> coordinates must be greater than "
-        "the <lower_left> coordinates on a tally mesh.");
-    }
-
-    // Set width
-    width_ = xt::eval((upper_right_ - lower_left_) / shape_);
-  } else {
-    fatal_error("Must specify either <upper_right> and <width> on a mesh.");
-  }
-
-  // Set volume fraction
-  volume_frac_ = 1.0/xt::prod(shape_)();
-}
-
-int RegularMesh::get_index_in_direction(double r, int i) const
-{
-  return std::ceil((r - lower_left_[i]) / width_[i]);
-}
-
-void RegularMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
-                               std::vector<double>& lengths) const
+void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
 {
   // ========================================================================
   // Determine where the track intersects the mesh and if it intersects at all.
@@ -565,10 +473,10 @@ void RegularMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
       if (std::fabs(u[k]) < FP_PRECISION) {
         d[k] = INFTY;
       } else if (u[k] > 0) {
-        double xyz_cross = lower_left_[k] + ijk0[k] * width_[k];
+        double xyz_cross = positive_grid_boundary(ijk0.data(), k);
         d[k] = (xyz_cross - r0[k]) / u[k];
       } else {
-        double xyz_cross = lower_left_[k] + (ijk0[k] - 1) * width_[k];
+        double xyz_cross = negative_grid_boundary(ijk0.data(), k);
         d[k] = (xyz_cross - r0[k]) / u[k];
       }
     }
@@ -602,6 +510,111 @@ void RegularMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
   }
 }
 
+
+//==============================================================================
+// RegularMesh implementation
+//==============================================================================
+
+RegularMesh::RegularMesh(pugi::xml_node node)
+  : StructuredMesh {node}
+{
+  // Determine number of dimensions for mesh
+  if (check_for_node(node, "dimension")) {
+    shape_ = get_node_xarray<int>(node, "dimension");
+    int n = n_dimension_ = shape_.size();
+    if (n != 1 && n != 2 && n != 3) {
+      fatal_error("Mesh must be one, two, or three dimensions.");
+    }
+
+    // Check that dimensions are all greater than zero
+    if (xt::any(shape_ <= 0)) {
+      fatal_error("All entries on the <dimension> element for a tally "
+        "mesh must be positive.");
+    }
+  }
+
+  // Check for lower-left coordinates
+  if (check_for_node(node, "lower_left")) {
+    // Read mesh lower-left corner location
+    lower_left_ = get_node_xarray<double>(node, "lower_left");
+  } else {
+    fatal_error("Must specify <lower_left> on a mesh.");
+  }
+
+  if (check_for_node(node, "width")) {
+    // Make sure both upper-right or width were specified
+    if (check_for_node(node, "upper_right")) {
+      fatal_error("Cannot specify both <upper_right> and <width> on a mesh.");
+    }
+
+    width_ = get_node_xarray<double>(node, "width");
+
+    // Check to ensure width has same dimensions
+    auto n = width_.size();
+    if (n != lower_left_.size()) {
+      fatal_error("Number of entries on <width> must be the same as "
+        "the number of entries on <lower_left>.");
+    }
+
+    // Check for negative widths
+    if (xt::any(width_ < 0.0)) {
+      fatal_error("Cannot have a negative <width> on a tally mesh.");
+    }
+
+    // Set width and upper right coordinate
+    upper_right_ = xt::eval(lower_left_ + shape_ * width_);
+
+  } else if (check_for_node(node, "upper_right")) {
+    upper_right_ = get_node_xarray<double>(node, "upper_right");
+
+    // Check to ensure width has same dimensions
+    auto n = upper_right_.size();
+    if (n != lower_left_.size()) {
+      fatal_error("Number of entries on <upper_right> must be the "
+        "same as the number of entries on <lower_left>.");
+    }
+
+    // Check that upper-right is above lower-left
+    if (xt::any(upper_right_ < lower_left_)) {
+      fatal_error("The <upper_right> coordinates must be greater than "
+        "the <lower_left> coordinates on a tally mesh.");
+    }
+
+    // Set width
+    if (shape_.size() > 0) {
+      width_ = xt::eval((upper_right_ - lower_left_) / shape_);
+    }
+  } else {
+    fatal_error("Must specify either <upper_right> and <width> on a mesh.");
+  }
+
+  // Make sure lower_left and dimension match
+  if (shape_.size() > 0) {
+    if (shape_.size() != lower_left_.size()) {
+      fatal_error("Number of entries on <lower_left> must be the same "
+        "as the number of entries on <dimension>.");
+    }
+
+    // Set volume fraction
+    volume_frac_ = 1.0/xt::prod(shape_)();
+  }
+}
+
+int RegularMesh::get_index_in_direction(double r, int i) const
+{
+  return std::ceil((r - lower_left_[i]) / width_[i]);
+}
+
+double RegularMesh::positive_grid_boundary(int* ijk, int i) const
+{
+  return lower_left_[i] + ijk[i] * width_[i];
+}
+
+double RegularMesh::negative_grid_boundary(int* ijk, int i) const
+{
+  return lower_left_[i] + (ijk[i] - 1) * width_[i];
+}
+
 void RegularMesh::surface_bins_crossed(const Particle& p,
                                        std::vector<int>& bins) const
 {
@@ -640,9 +653,9 @@ void RegularMesh::surface_bins_crossed(const Particle& p,
   Position xyz_cross;
   for (int i = 0; i < n; ++i) {
     if (u[i] > 0.0) {
-      xyz_cross[i] = lower_left_[i] + ijk0[i] * width_[i];
+      xyz_cross[i] = positive_grid_boundary(ijk0.data(), i);
     } else {
-      xyz_cross[i] = lower_left_[i] + (ijk0[i] - 1) * width_[i];
+      xyz_cross[i] = negative_grid_boundary(ijk0.data(), i);
     }
   }
 
@@ -887,109 +900,14 @@ RectilinearMesh::RectilinearMesh(pugi::xml_node node)
   upper_right_ = {grid_[0].back(), grid_[1].back(), grid_[2].back()};
 }
 
-void RectilinearMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
-                                   std::vector<double>& lengths) const
+double RectilinearMesh::positive_grid_boundary(int* ijk, int i) const
 {
-  // ========================================================================
-  // Determine where the track intersects the mesh and if it intersects at all.
+  return grid_[i][ijk[i]];
+}
 
-  // Copy the starting and ending coordinates of the particle.
-  Position last_r {p.r_last_};
-  Position r {p.r()};
-  Direction u {p.u()};
-
-  // Compute the length of the entire track.
-  double total_distance = (r - last_r).norm();
-
-  // While determining if this track intersects the mesh, offset the starting
-  // and ending coords by a bit.  This avoid finite-precision errors that can
-  // occur when the mesh surfaces coincide with lattice or geometric surfaces.
-  Position r0 = last_r + TINY_BIT*u;
-  Position r1 = r - TINY_BIT*u;
-
-  // Determine the mesh indices for the starting and ending coords.
-  int ijk0[3], ijk1[3];
-  bool start_in_mesh;
-  get_indices(r0, ijk0, &start_in_mesh);
-  bool end_in_mesh;
-  get_indices(r1, ijk1, &end_in_mesh);
-
-  // Reset coordinates and check for a mesh intersection if necessary.
-  if (start_in_mesh) {
-    // The initial coords lie in the mesh, use those coords for tallying.
-    r0 = last_r;
-  } else {
-    // The initial coords do not lie in the mesh.  Check to see if the particle
-    // eventually intersects the mesh and compute the relevant coords and
-    // indices.
-    if (!intersects(r0, r1, ijk0)) return;
-  }
-  r1 = r;
-
-  // The TINY_BIT offsets above mean that the preceding logic cannot always find
-  // the correct ijk0 and ijk1 indices. For tracks shorter than 2*TINY_BIT, just
-  // assume the track lies in only one mesh bin. These tracks are very short so
-  // any error caused by this assumption will be small. It is important that
-  // ijk0 values are used rather than ijk1 because the previous logic guarantees
-  // ijk0 is a valid mesh bin.
-  if (total_distance < 2*TINY_BIT) {
-    for (int i = 0; i < 3; ++i) ijk1[i] = ijk0[i];
-  }
-
-  // ========================================================================
-  // Find which mesh cells are traversed and the length of each traversal.
-
-  while (true) {
-    if (std::equal(ijk0, ijk0+3, ijk1)) {
-      // The track ends in this cell.  Use the particle end location rather
-      // than the mesh surface and stop iterating.
-      double distance = (r1 - r0).norm();
-      bins.push_back(get_bin_from_indices(ijk0));
-      lengths.push_back(distance / total_distance);
-      break;
-    }
-
-    // The track exits this cell.  Determine the distance to each mesh surface.
-    double d[3];
-    for (int k = 0; k < 3; ++k) {
-      if (std::fabs(u[k]) < FP_PRECISION) {
-        d[k] = INFTY;
-      } else if (u[k] > 0) {
-        double xyz_cross = grid_[k][ijk0[k]];
-        d[k] = (xyz_cross - r0[k]) / u[k];
-      } else {
-        double xyz_cross = grid_[k][ijk0[k] - 1];
-        d[k] = (xyz_cross - r0[k]) / u[k];
-      }
-    }
-
-    // Pick the closest mesh surface and append this traversal to the output.
-    auto j = std::min_element(d, d+3) - d;
-    double distance = d[j];
-    bins.push_back(get_bin_from_indices(ijk0));
-    lengths.push_back(distance / total_distance);
-
-    // Translate to the oncoming mesh surface.
-    r0 += distance * u;
-
-    // Increment the indices into the next mesh cell.
-    if (u[j] > 0.0) {
-      ++ijk0[j];
-    } else {
-      --ijk0[j];
-    }
-
-    // If the next indices are invalid, then the track has left the mesh and
-    // we are done.
-    bool in_mesh = true;
-    for (int i = 0; i < 3; ++i) {
-      if (ijk0[i] < 1 || ijk0[i] > shape_[i]) {
-        in_mesh = false;
-        break;
-      }
-    }
-    if (!in_mesh) break;
-  }
+double RectilinearMesh::negative_grid_boundary(int* ijk, int i) const
+{
+  return grid_[i][ijk[i] - 1];
 }
 
 void RectilinearMesh::surface_bins_crossed(const Particle& p,
@@ -1086,9 +1004,9 @@ void RectilinearMesh::surface_bins_crossed(const Particle& p,
   Position xyz_cross;
   for (int i = 0; i < 3; ++i) {
     if (u[i] > 0.0) {
-      xyz_cross[i] = grid_[i][ijk0[i]];
+      xyz_cross[i] = positive_grid_boundary(ijk0, i);
     } else {
-      xyz_cross[i] = grid_[i][ijk0[i] - 1];
+      xyz_cross[i] = negative_grid_boundary(ijk0, i);
     }
   }
 

From 87c22dd4d498b6c28621c690b12d75da37fd1665 Mon Sep 17 00:00:00 2001
From: aprilnovak <novak@berkeley.edu>
Date: Mon, 19 Oct 2020 16:19:28 -0500
Subject: [PATCH 8/8] Use arrays for ijk0, ijk1, and d in bins_crossed for
 slightly better performance. Refs #1695

---
 src/mesh.cpp | 30 ++++++++++++++++++------------
 1 file changed, 18 insertions(+), 12 deletions(-)

diff --git a/src/mesh.cpp b/src/mesh.cpp
index d3b671263a..53279cc57a 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -406,6 +406,7 @@ bool StructuredMesh::intersects_3d(Position& r0, Position r1, int* ijk) const
 }
 
 void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
+  std::vector<double>& lengths) const
 {
   // ========================================================================
   // Determine where the track intersects the mesh and if it intersects at all.
@@ -424,13 +425,18 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
   Position r0 = last_r + TINY_BIT*u;
   Position r1 = r - TINY_BIT*u;
 
-  // Determine the mesh indices for the starting and ending coords.
+  // Determine the mesh indices for the starting and ending coords. Here, we
+  // use arrays for ijk0 and ijk1 instead of std::vector because we obtain a
+  // small performance improvement by forcing this data to live on the stack,
+  // rather than on the heap. We know the maximum length is 3, and by
+  // ensuring that all loops are only indexed up to n_dimension, we will not
+  // access any non-initialized values. The same concept is used throughout.
   int n = n_dimension_;
-  std::vector<int> ijk0(n), ijk1(n);
+  int ijk0[3], ijk1[3];
   bool start_in_mesh;
-  get_indices(r0, ijk0.data(), &start_in_mesh);
+  get_indices(r0, ijk0, &start_in_mesh);
   bool end_in_mesh;
-  get_indices(r1, ijk1.data(), &end_in_mesh);
+  get_indices(r1, ijk1, &end_in_mesh);
 
   // Reset coordinates and check for a mesh intersection if necessary.
   if (start_in_mesh) {
@@ -440,7 +446,7 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
     // The initial coords do not lie in the mesh.  Check to see if the particle
     // eventually intersects the mesh and compute the relevant coords and
     // indices.
-    if (!intersects(r0, r1, ijk0.data())) return;
+    if (!intersects(r0, r1, ijk0)) return;
   }
   r1 = r;
 
@@ -458,33 +464,33 @@ void StructuredMesh::bins_crossed(const Particle& p, std::vector<int>& bins,
   // Find which mesh cells are traversed and the length of each traversal.
 
   while (true) {
-    if (ijk0 == ijk1) {
+    if (std::equal(ijk0, ijk0 + n, ijk1)) {
       // The track ends in this cell.  Use the particle end location rather
       // than the mesh surface and stop iterating.
       double distance = (r1 - r0).norm();
-      bins.push_back(get_bin_from_indices(ijk0.data()));
+      bins.push_back(get_bin_from_indices(ijk0));
       lengths.push_back(distance / total_distance);
       break;
     }
 
     // The track exits this cell.  Determine the distance to each mesh surface.
-    std::vector<double> d(n);
+    double d[3];
     for (int k = 0; k < n; ++k) {
       if (std::fabs(u[k]) < FP_PRECISION) {
         d[k] = INFTY;
       } else if (u[k] > 0) {
-        double xyz_cross = positive_grid_boundary(ijk0.data(), k);
+        double xyz_cross = positive_grid_boundary(ijk0, k);
         d[k] = (xyz_cross - r0[k]) / u[k];
       } else {
-        double xyz_cross = negative_grid_boundary(ijk0.data(), k);
+        double xyz_cross = negative_grid_boundary(ijk0, k);
         d[k] = (xyz_cross - r0[k]) / u[k];
       }
     }
 
     // Pick the closest mesh surface and append this traversal to the output.
-    auto j = std::min_element(d.begin(), d.end()) - d.begin();
+    auto j = std::min_element(d, d + n) - d;
     double distance = d[j];
-    bins.push_back(get_bin_from_indices(ijk0.data()));
+    bins.push_back(get_bin_from_indices(ijk0));
     lengths.push_back(distance / total_distance);
 
     // Translate to the oncoming mesh surface.