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