Fix get_index_in_direction for regular meshes

include/openmc/mesh.h

@@ -374,6 +374,9 @@ class StructuredMesh : public Mesh {
   //!
   //! \param[in] r Coordinate to get index for
   //! \param[in] i Direction index
+  //! \return Mesh index in [0, shape[i] + 1]. The external boundaries are
+  //! included in the mesh, and interior boundaries belong to the lower-index
+  //! mesh cell.
   virtual int get_index_in_direction(double r, int i) const = 0;
 
   //! Get the coordinate for the mesh grid boundary in the positive direction

src/mesh.cpp

@@ -1472,6 +1472,11 @@ RegularMesh::RegularMesh(hid_t group) : StructuredMesh {group}
 
 int RegularMesh::get_index_in_direction(double r, int i) const
 {
+  if (r <= lower_left_[i])
+    return r == lower_left_[i] ? 1 : 0;
+  if (r >= upper_right_[i])
+    return r == upper_right_[i] ? shape_[i] : shape_[i] + 1;
+
   return std::ceil((r - lower_left_[i]) / width_[i]);
 }
 

tests/cpp_unit_tests/test_mesh.cpp

@@ -2,6 +2,7 @@
 #include <iostream>
 #include <string>
 
+#include <catch2/catch_approx.hpp>
 #include <catch2/catch_test_macros.hpp>
 #include <pugixml.hpp>
 
@@ -255,3 +256,116 @@ TEST_CASE("Test multiple meshes HDF5 roundtrip - spherical")
   REQUIRE(regular_mesh_hdf5->lower_left() == regular_mesh_xml->lower_left());
   REQUIRE(regular_mesh_hdf5->upper_right() == regular_mesh_xml->upper_right());
 }
+
+TEST_CASE("Test get_index_in_direction - regular")
+{
+  // The XML data as a string
+  std::string xml_string = R"(
+        <mesh id="1">
+            <dimension>2 1 1</dimension>
+            <lower_left>-2 -2 -2</lower_left>
+            <upper_right>2 2 2</upper_right>
+       </mesh>
+    )";
+
+  // Create a pugixml document object
+  pugi::xml_document doc;
+
+  // Load the XML from the string
+  pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
+
+  pugi::xml_node root = doc.child("mesh");
+
+  auto mesh = RegularMesh(root);
+
+  REQUIRE(mesh.get_index_in_direction(-4.1, 0) == 0);
+  REQUIRE(mesh.get_index_in_direction(-4.0, 0) == 0);
+  REQUIRE(mesh.get_index_in_direction(-3.9, 0) == 0);
+
+  REQUIRE(mesh.get_index_in_direction(-2.1, 0) == 0);
+  REQUIRE(mesh.get_index_in_direction(-2.0, 0) == 1); // lower left
+  REQUIRE(mesh.get_index_in_direction(-1.9, 0) == 1);
+
+  REQUIRE(mesh.get_index_in_direction(-0.1, 0) == 1);
+  REQUIRE(mesh.get_index_in_direction(0.0, 0) == 1);
+  REQUIRE(mesh.get_index_in_direction(0.1, 0) == 2);
+
+  REQUIRE(mesh.get_index_in_direction(1.9, 0) == 2);
+  REQUIRE(mesh.get_index_in_direction(2.0, 0) == 2); // upper right
+  REQUIRE(mesh.get_index_in_direction(2.1, 0) == 3);
+
+  REQUIRE(mesh.get_index_in_direction(3.9, 0) == 3);
+  REQUIRE(mesh.get_index_in_direction(4.0, 0) == 3);
+  REQUIRE(mesh.get_index_in_direction(4.1, 0) == 3);
+}
+
+TEST_CASE("Test get_index_in_direction - rectilinear")
+{
+  // The XML data as a string
+  std::string xml_string = R"(
+        <mesh id="1" type="rectilinear">
+            <x_grid>-1.0 0.0 2.0</x_grid>
+            <y_grid>-1.0 1.0</y_grid>
+            <z_grid>-1.0 1.0</z_grid>
+        </mesh>
+    )";
+
+  // Create a pugixml document object
+  pugi::xml_document doc;
+
+  // Load the XML from the string
+  pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
+
+  pugi::xml_node root = doc.child("mesh");
+
+  auto mesh = RectilinearMesh(root);
+
+  REQUIRE(mesh.get_index_in_direction(-5.0, 0) == 0);
+
+  REQUIRE(mesh.get_index_in_direction(-1.1, 0) == 0);
+  REQUIRE(mesh.get_index_in_direction(-1.0, 0) == 1); // lower left
+  REQUIRE(mesh.get_index_in_direction(-0.9, 0) == 1);
+
+  REQUIRE(mesh.get_index_in_direction(1.9, 0) == 2);
+  REQUIRE(mesh.get_index_in_direction(2.0, 0) == 2); // upper right
+  REQUIRE(mesh.get_index_in_direction(2.1, 0) == 3);
+
+  REQUIRE(mesh.get_index_in_direction(5.0, 0) == 3);
+}
+
+TEST_CASE("Test regular mesh ray tracing from outside")
+{
+  std::string xml_string = R"(
+        <mesh id="1">
+            <dimension>2 1 1</dimension>
+            <lower_left>-2 -2 -2</lower_left>
+            <upper_right>2 2 2</upper_right>
+       </mesh>
+    )";
+
+  pugi::xml_document doc;
+  pugi::xml_parse_result result = doc.load_string(xml_string.c_str());
+  auto mesh = RegularMesh(doc.child("mesh"));
+
+  vector<int> bins;
+  vector<double> lengths;
+
+  mesh.bins_crossed(Position {-6.0, 0.0, 0.0}, Position {6.0, 0.0, 0.0},
+    Direction {1.0, 0.0, 0.0}, bins, lengths);
+
+  REQUIRE(bins == vector<int> {0, 1});
+  REQUIRE(lengths.size() == 2);
+  REQUIRE(lengths[0] == Catch::Approx(1.0 / 6.0));
+  REQUIRE(lengths[1] == Catch::Approx(1.0 / 6.0));
+
+  bins.clear();
+  lengths.clear();
+
+  mesh.bins_crossed(Position {6.0, 0.0, 0.0}, Position {-6.0, 0.0, 0.0},
+    Direction {-1.0, 0.0, 0.0}, bins, lengths);
+
+  REQUIRE(bins == vector<int> {1, 0});
+  REQUIRE(lengths.size() == 2);
+  REQUIRE(lengths[0] == Catch::Approx(1.0 / 6.0));
+  REQUIRE(lengths[1] == Catch::Approx(1.0 / 6.0));
+}

tests/regression_tests/weightwindows/survival_biasing/local/results_true.dat

@@ -1 +1 @@
-386e507008ed3c72c6e1a101aafc01cfaaff2c0b10555558d1eab6332441f7b17264b55002d721151bac2f3e7c1a8aac4c5ed243f5270533d171850f985206ed
+0d7b17d4e364bda2be21feb2e5b2aae4be69fc8ed634c4f45062e8efc61b7bd24dcf448837692fd603afe3756501fe8821d134f2d6289179618f85178ddb8793

tests/regression_tests/weightwindows/survival_biasing/shared/results_true.dat

@@ -1 +1 @@
-d17a437262d3316985fba4b48e21a7fcd5f32ac2d96ef0e66849f567deaeadc1e0f18d5d0bf5e9cab4246dbf823e7f43f249a1f67e928b27ae8c70b89f3cefbf
+2f1c8efbf6ab0b0c7319fadd7ee8f27b545c1ebab94ebe15cf1133b54227487542966b6d798970c8b37bc19e7041e687687d5285803241bb83f4c8dc40a4d276