Prevent segfault in distance to boundary calculation

src/geometry.cpp

@@ -419,13 +419,14 @@ BoundaryInfo distance_to_boundary(Particle& p)
     double& d = info.distance;
     if (d_surf < d_lat - FP_COINCIDENT) {
       if (d == INFINITY || (d - d_surf) / d >= FP_REL_PRECISION) {
+        // Update closest distance
         d = d_surf;
 
         // If the cell is not simple, it is possible that both the negative and
         // positive half-space were given in the region specification. Thus, we
         // have to explicitly check which half-space the particle would be
         // traveling into if the surface is crossed
-        if (c.is_simple()) {
+        if (c.is_simple() || d == INFTY) {
           info.surface_index = level_surf_cross;
         } else {
           Position r_hit = r + d_surf * u;

tests/unit_tests/test_no_visible_boundary.py

@@ -0,0 +1,30 @@
+import openmc
+
+
+def test_no_visible_boundary(run_in_tmpdir):
+    copper = openmc.Material()
+    copper.add_nuclide('Cu63', 1.0)
+    copper.set_density('g/cm3', 0.3)
+    air = openmc.Material()
+    air.add_nuclide('N14', 1.0)
+    air.set_density('g/cm3', 0.0012)
+
+    # Create a simple model of a neutron source directly impinging on a thin
+    # disc of copper. Neutrons leaving the back of the disc see no surfaces in
+    # front of them.
+    disc = openmc.model.RightCircularCylinder((0., 0., 1.), 0.1, 1.2)
+    box = openmc.rectangular_prism(width=10, height=10, boundary_type='vacuum')
+    c1 = openmc.Cell(fill=copper, region=-disc)
+    c2 = openmc.Cell(fill=air, region=+disc & box)
+    model = openmc.Model()
+    model.geometry = openmc.Geometry([c1, c2])
+    model.settings.run_mode = 'fixed source'
+    model.settings.particles = 1000
+    model.settings.batches = 5
+    model.settings.source = openmc.IndependentSource(
+        space=openmc.stats.Point(),
+        angle=openmc.stats.Monodirectional((0., 0., 1.))
+    )
+
+    # Run model to ensure it doesn't segfault
+    model.run()