diff --git a/src/mesh.cpp b/src/mesh.cpp index acd80c7ba..9e27fc4e6 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -524,6 +524,14 @@ void RegularMesh::bins_crossed(const Particle* p, std::vector& bins, } 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. + if (total_distance < 2*TINY_BIT) { + for (int i = 0; i < n; ++i) ijk0[i] = ijk1[i]; + } + // ======================================================================== // Find which mesh cells are traversed and the length of each traversal. @@ -899,6 +907,14 @@ void RectilinearMesh::bins_crossed(const Particle* p, std::vector& bins, } 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. + if (total_distance < 2*TINY_BIT) { + for (int i = 0; i < 3; ++i) ijk0[i] = ijk1[i]; + } + // ======================================================================== // Find which mesh cells are traversed and the length of each traversal.