Refactor Ray class into its own file

CMakeLists.txt

@@ -407,6 +407,7 @@ list(APPEND libopenmc_SOURCES
   src/random_ray/linear_source_domain.cpp
   src/random_ray/moment_matrix.cpp
   src/random_ray/source_region.cpp
+  src/ray.cpp
   src/reaction.cpp
   src/reaction_product.cpp
   src/scattdata.cpp

include/openmc/plot.h

@@ -17,6 +17,7 @@
 #include "openmc/particle.h"
 #include "openmc/position.h"
 #include "openmc/random_lcg.h"
+#include "openmc/ray.h"
 #include "openmc/xml_interface.h"
 
 namespace openmc {
@@ -497,47 +498,6 @@ class SolidRayTracePlot : public RayTracePlot {
   Position light_location_;
 };
 
-// Base class that implements ray tracing logic, not necessarily through
-// defined regions of the geometry but also outside of it.
-class Ray : public GeometryState {
-
-public:
-  // Initialize from location and direction
-  Ray(Position r, Direction u) { init_from_r_u(r, u); }
-
-  // Initialize from known geometry state
-  Ray(const GeometryState& p) : GeometryState(p) {}
-
-  // Called at every surface intersection within the model
-  virtual void on_intersection() = 0;
-
-  /*
-   * Traces the ray through the geometry, calling on_intersection
-   * at every surface boundary.
-   */
-  void trace();
-
-  // Stops the ray and exits tracing when called from on_intersection
-  void stop() { stop_ = true; }
-
-  // Sets the dist_ variable
-  void compute_distance();
-
-protected:
-  // Records how far the ray has traveled
-  double traversal_distance_ {0.0};
-
-private:
-  // Max intersections before we assume ray tracing is caught in an infinite
-  // loop:
-  static const int MAX_INTERSECTIONS = 1000000;
-
-  bool hit_something_ {false};
-  bool stop_ {false};
-
-  unsigned event_counter_ {0};
-};
-
 class ProjectionRay : public Ray {
 public:
   ProjectionRay(Position r, Direction u, const WireframeRayTracePlot& plot,

include/openmc/ray.h

@@ -0,0 +1,50 @@
+#ifndef OPENMC_RAY_H
+#define OPENMC_RAY_H
+
+#include "openmc/particle_data.h"
+#include "openmc/position.h"
+
+namespace openmc {
+
+// Base class that implements ray tracing logic, not necessarily through
+// defined regions of the geometry but also outside of it.
+class Ray : public GeometryState {
+
+public:
+  // Initialize from location and direction
+  Ray(Position r, Direction u) { init_from_r_u(r, u); }
+
+  // Initialize from known geometry state
+  Ray(const GeometryState& p) : GeometryState(p) {}
+
+  // Called at every surface intersection within the model
+  virtual void on_intersection() = 0;
+
+  /*
+   * Traces the ray through the geometry, calling on_intersection
+   * at every surface boundary.
+   */
+  void trace();
+
+  // Stops the ray and exits tracing when called from on_intersection
+  void stop() { stop_ = true; }
+
+  // Sets the dist_ variable
+  void compute_distance();
+
+protected:
+  // Records how far the ray has traveled
+  double traversal_distance_ {0.0};
+
+private:
+  // Max intersections before we assume ray tracing is caught in an infinite
+  // loop:
+  static const int MAX_INTERSECTIONS = 1000000;
+
+  bool stop_ {false};
+
+  unsigned event_counter_ {0};
+};
+
+} // namespace openmc
+#endif // OPENMC_RAY_H

src/plot.cpp

@@ -1649,165 +1649,6 @@ void SolidRayTracePlot::set_diffuse_fraction(pugi::xml_node node)
   }
 }
 
-void Ray::compute_distance()
-{
-  boundary() = distance_to_boundary(*this);
-}
-
-void Ray::trace()
-{
-  // To trace the ray from its origin all the way through the model, we have
-  // to proceed in two phases. In the first, the ray may or may not be found
-  // inside the model. If the ray is already in the model, phase one can be
-  // skipped. Otherwise, the ray has to be advanced to the boundary of the
-  // model where all the cells are defined. Importantly, this is assuming that
-  // the model is convex, which is a very reasonable assumption for any
-  // radiation transport model.
-  //
-  // After phase one is done, we can starting tracing from cell to cell within
-  // the model. This step can use neighbor lists to accelerate the ray tracing.
-
-  bool inside_cell;
-  // Check for location if the particle is already known
-  if (lowest_coord().cell() == C_NONE) {
-    // The geometry position of the particle is either unknown or outside of the
-    // edge of the model.
-    if (lowest_coord().universe() == C_NONE) {
-      // Attempt to initialize the particle. We may have to
-      // enter a loop to move it up to the edge of the model.
-      inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
-    } else {
-      // It has been already calculated that the current position is outside of
-      // the edge of the model.
-      inside_cell = false;
-    }
-  } else {
-    // Availability of the cell means that the particle is located inside the
-    // edge.
-    inside_cell = true;
-  }
-
-  // Advance to the boundary of the model
-  while (!inside_cell) {
-    advance_to_boundary_from_void();
-    inside_cell = exhaustive_find_cell(*this, settings::verbosity >= 10);
-
-    // If true this means no surface was intersected. See cell.cpp and search
-    // for numeric_limits to see where we return it.
-    if (surface() == std::numeric_limits<int>::max()) {
-      warning(fmt::format("Lost a ray, r = {}, u = {}", r(), u()));
-      return;
-    }
-
-    // Exit this loop and enter into cell-to-cell ray tracing (which uses
-    // neighbor lists)
-    if (inside_cell)
-      break;
-
-    // if there is no intersection with the model, we're done
-    if (boundary().surface() == SURFACE_NONE)
-      return;
-
-    event_counter_++;
-    if (event_counter_ > MAX_INTERSECTIONS) {
-      warning("Likely infinite loop in ray traced plot");
-      return;
-    }
-  }
-
-  // Call the specialized logic for this type of ray. This is for the
-  // intersection for the first intersection if we had one.
-  if (boundary().surface() != SURFACE_NONE) {
-    // set the geometry state's surface attribute to be used for
-    // surface normal computation
-    surface() = boundary().surface();
-    on_intersection();
-    if (stop_)
-      return;
-  }
-
-  // reset surface attribute to zero after the first intersection so that it
-  // doesn't perturb surface crossing logic from here on out
-  surface() = 0;
-
-  // This is the ray tracing loop within the model. It exits after exiting
-  // the model, which is equivalent to assuming that the model is convex.
-  // It would be nice to factor out the on_intersection at the end of this
-  // loop and then do "while (inside_cell)", but we can't guarantee it's
-  // on a surface in that case. There might be some other way to set it
-  // up that is perhaps a little more elegant, but this is what works just
-  // fine.
-  while (true) {
-
-    compute_distance();
-
-    // There are no more intersections to process
-    // if we hit the edge of the model, so stop
-    // the particle in that case. Also, just exit
-    // if a negative distance was somehow computed.
-    if (boundary().distance() == INFTY || boundary().distance() == INFINITY ||
-        boundary().distance() < 0) {
-      return;
-    }
-
-    // See below comment where call_on_intersection is checked in an
-    // if statement for an explanation of this.
-    bool call_on_intersection {true};
-    if (boundary().distance() < 10 * TINY_BIT) {
-      call_on_intersection = false;
-    }
-
-    // DAGMC surfaces expect us to go a little bit further than the advance
-    // distance to properly check cell inclusion.
-    boundary().distance() += TINY_BIT;
-
-    // Advance particle, prepare for next intersection
-    for (int lev = 0; lev < n_coord(); ++lev) {
-      coord(lev).r() += boundary().distance() * coord(lev).u();
-    }
-    surface() = boundary().surface();
-    // Initialize last cells from the current cell, because the cell() variable
-    // does not contain the data for the case of a single-segment ray
-    for (int j = 0; j < n_coord(); ++j) {
-      cell_last(j) = coord(j).cell();
-    }
-    n_coord_last() = n_coord();
-    n_coord() = boundary().coord_level();
-    if (boundary().lattice_translation()[0] != 0 ||
-        boundary().lattice_translation()[1] != 0 ||
-        boundary().lattice_translation()[2] != 0) {
-      cross_lattice(*this, boundary(), settings::verbosity >= 10);
-    }
-
-    // Record how far the ray has traveled
-    traversal_distance_ += boundary().distance();
-    inside_cell = neighbor_list_find_cell(*this, settings::verbosity >= 10);
-
-    // Call the specialized logic for this type of ray. Note that we do not
-    // call this if the advance distance is very small. Unfortunately, it seems
-    // darn near impossible to get the particle advanced to the model boundary
-    // and through it without sometimes accidentally calling on_intersection
-    // twice. This incorrectly shades the region as occluded when it might not
-    // actually be. By screening out intersection distances smaller than a
-    // threshold 10x larger than the scoot distance used to advance up to the
-    // model boundary, we can avoid that situation.
-    if (call_on_intersection) {
-      on_intersection();
-      if (stop_)
-        return;
-    }
-
-    if (!inside_cell)
-      return;
-
-    event_counter_++;
-    if (event_counter_ > MAX_INTERSECTIONS) {
-      warning("Likely infinite loop in ray traced plot");
-      return;
-    }
-  }
-}
-
 void ProjectionRay::on_intersection()
 {
   // This records a tuple with the following info