Add constant for invalid surface tokens. (#3260)

Co-authored-by: Paul Romano <[email protected]>

include/openmc/cell.h

@@ -29,7 +29,6 @@ namespace openmc {
 
 enum class Fill { MATERIAL, UNIVERSE, LATTICE };
 
-// TODO: Convert to enum
 constexpr int32_t OP_LEFT_PAREN {std::numeric_limits<int32_t>::max()};
 constexpr int32_t OP_RIGHT_PAREN {std::numeric_limits<int32_t>::max() - 1};
 constexpr int32_t OP_COMPLEMENT {std::numeric_limits<int32_t>::max() - 2};

include/openmc/constants.h

@@ -350,6 +350,10 @@ enum class RandomRaySourceShape { FLAT, LINEAR, LINEAR_XY };
 
 enum class GeometryType { CSG, DAG };
 
+// a surface token cannot be zero due to the unsigned nature of zero for integer
+// representations. This value represents no surface.
+constexpr int32_t SURFACE_NONE {0};
+
 } // namespace openmc
 
 #endif // OPENMC_CONSTANTS_H

include/openmc/particle_data.h

@@ -185,10 +185,14 @@ struct CacheDataMG {
 
 struct BoundaryInfo {
   double distance {INFINITY}; //!< distance to nearest boundary
-  int surface_index {0}; //!< if boundary is surface, index in surfaces vector
-  int coord_level;       //!< coordinate level after crossing boundary
+  int surface {
+    SURFACE_NONE}; //!< surface token, non-zero if boundary is surface
+  int coord_level; //!< coordinate level after crossing boundary
   array<int, 3>
     lattice_translation {}; //!< which way lattice indices will change
+
+  // TODO: off-by-one
+  int surface_index() const { return std::abs(surface) - 1; }
 };
 
 /*
@@ -226,7 +230,7 @@ class GeometryState {
   void init_from_r_u(Position r_a, Direction u_a)
   {
     clear();
-    surface() = 0;
+    surface() = SURFACE_NONE;
     material() = C_NONE;
     r() = r_a;
     u() = u_a;
@@ -296,10 +300,17 @@ class GeometryState {
   Direction& u_local() { return coord_[n_coord_ - 1].u; }
   const Direction& u_local() const { return coord_[n_coord_ - 1].u; }
 
-  // Surface that the particle is on
+  // Surface token for the surface that the particle is currently on
   int& surface() { return surface_; }
   const int& surface() const { return surface_; }
 
+  // Surface index based on the current value of the surface_ attribute
+  int surface_index() const
+  {
+    // TODO: off-by-one
+    return std::abs(surface_) - 1;
+  }
+
   // Boundary information
   BoundaryInfo& boundary() { return boundary_; }
 
@@ -337,7 +348,8 @@ class GeometryState {
   Position r_last_;         //!< previous coordinates
   Direction u_last_;        //!< previous direction coordinates
 
-  int surface_ {0}; //!< index for surface particle is on
+  int surface_ {
+    SURFACE_NONE}; //!< surface token for surface the particle is currently on
 
   BoundaryInfo boundary_; //!< Info about the next intersection
 

src/boundary_condition.cpp

@@ -129,8 +129,7 @@ TranslationalPeriodicBC::TranslationalPeriodicBC(int i_surf, int j_surf)
 void TranslationalPeriodicBC::handle_particle(
   Particle& p, const Surface& surf) const
 {
-  // TODO: off-by-one on surface indices throughout this function.
-  int i_particle_surf = std::abs(p.surface()) - 1;
+  int i_particle_surf = p.surface_index();
 
   // Figure out which of the two BC surfaces were struck then find the
   // particle's new location and surface.
@@ -255,8 +254,7 @@ RotationalPeriodicBC::RotationalPeriodicBC(int i_surf, int j_surf)
 void RotationalPeriodicBC::handle_particle(
   Particle& p, const Surface& surf) const
 {
-  // TODO: off-by-one on surface indices throughout this function.
-  int i_particle_surf = std::abs(p.surface()) - 1;
+  int i_particle_surf = p.surface_index();
 
   // Figure out which of the two BC surfaces were struck to figure out if a
   // forward or backward rotation is required.  Specify the other surface as

src/dagmc.cpp

@@ -847,7 +847,7 @@ void check_dagmc_root_univ()
 
 int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ)
 {
-  auto surfp = dynamic_cast<DAGSurface*>(model::surfaces[surf - 1].get());
+  auto surfp = dynamic_cast<DAGSurface*>(model::surfaces[surf].get());
   auto cellp = dynamic_cast<DAGCell*>(model::cells[curr_cell].get());
   auto univp = static_cast<DAGUniverse*>(model::universes[univ].get());
 

src/geometry.cpp

@@ -421,15 +421,15 @@ BoundaryInfo distance_to_boundary(GeometryState& p)
         // have to explicitly check which half-space the particle would be
         // traveling into if the surface is crossed
         if (c.is_simple() || d == INFTY) {
-          info.surface_index = level_surf_cross;
+          info.surface = level_surf_cross;
         } else {
           Position r_hit = r + d_surf * u;
           Surface& surf {*model::surfaces[std::abs(level_surf_cross) - 1]};
           Direction norm = surf.normal(r_hit);
           if (u.dot(norm) > 0) {
-            info.surface_index = std::abs(level_surf_cross);
+            info.surface = std::abs(level_surf_cross);
           } else {
-            info.surface_index = -std::abs(level_surf_cross);
+            info.surface = -std::abs(level_surf_cross);
           }
         }
 
@@ -441,7 +441,7 @@ BoundaryInfo distance_to_boundary(GeometryState& p)
     } else {
       if (d == INFINITY || (d - d_lat) / d >= FP_REL_PRECISION) {
         d = d_lat;
-        info.surface_index = 0;
+        info.surface = SURFACE_NONE;
         info.lattice_translation = level_lat_trans;
         info.coord_level = i + 1;
       }

src/output.cpp

@@ -200,9 +200,9 @@ void print_particle(Particle& p)
   }
 
   // Display miscellaneous info.
-  if (p.surface() != 0) {
+  if (p.surface() != SURFACE_NONE) {
     // Surfaces identifiers are >= 1, but indices are >= 0 so we need -1
-    const Surface& surf {*model::surfaces[std::abs(p.surface()) - 1]};
+    const Surface& surf {*model::surfaces[p.surface_index()]};
     fmt::print("  Surface = {}\n", (p.surface() > 0) ? surf.id_ : -surf.id_);
   }
   fmt::print("  Weight = {}\n", p.wgt());

src/particle.cpp

@@ -108,7 +108,7 @@ void Particle::from_source(const SourceSite* src)
 {
   // Reset some attributes
   clear();
-  surface() = 0;
+  surface() = SURFACE_NONE;
   cell_born() = C_NONE;
   material() = C_NONE;
   n_collision() = 0;
@@ -277,7 +277,7 @@ void Particle::event_cross_surface()
   n_coord_last() = n_coord();
 
   // Set surface that particle is on and adjust coordinate levels
-  surface() = boundary().surface_index;
+  surface() = boundary().surface;
   n_coord() = boundary().coord_level;
 
   if (boundary().lattice_translation[0] != 0 ||
@@ -291,7 +291,7 @@ void Particle::event_cross_surface()
   } else {
     // Particle crosses surface
     // TODO: off-by-one
-    const auto& surf {model::surfaces[std::abs(surface()) - 1].get()};
+    const auto& surf {model::surfaces[surface_index()].get()};
     // If BC, add particle to surface source before crossing surface
     if (surf->surf_source_ && surf->bc_) {
       add_surf_source_to_bank(*this, *surf);
@@ -328,7 +328,7 @@ void Particle::event_collide()
     score_surface_tally(*this, model::active_meshsurf_tallies);
 
   // Clear surface component
-  surface() = 0;
+  surface() = SURFACE_NONE;
 
   if (settings::run_CE) {
     collision(*this);
@@ -549,7 +549,7 @@ void Particle::cross_surface(const Surface& surf)
 #ifdef DAGMC
   // in DAGMC, we know what the next cell should be
   if (surf.geom_type() == GeometryType::DAG) {
-    int32_t i_cell = next_cell(std::abs(surface()), cell_last(n_coord() - 1),
+    int32_t i_cell = next_cell(surface_index(), cell_last(n_coord() - 1),
                        lowest_coord().universe) -
                      1;
     // save material and temp
@@ -587,7 +587,7 @@ void Particle::cross_surface(const Surface& surf)
     // the particle is really traveling tangent to a surface, if we move it
     // forward a tiny bit it should fix the problem.
 
-    surface() = 0;
+    surface() = SURFACE_NONE;
     n_coord() = 1;
     r() += TINY_BIT * u();
 

src/plot.cpp

@@ -1301,7 +1301,7 @@ void ProjectionPlot::create_output() const
           while (intersection_found) {
             bool inside_cell = false;
 
-            int32_t i_surface = std::abs(p.surface()) - 1;
+            int32_t i_surface = p.surface_index();
             if (i_surface > 0 &&
                 model::surfaces[i_surface]->geom_type() == GeometryType::DAG) {
 #ifdef DAGMC
@@ -1334,13 +1334,13 @@ void ProjectionPlot::create_output() const
               this_line_segments[tid][horiz].emplace_back(
                 color_by_ == PlotColorBy::mats ? p.material()
                                                : p.lowest_coord().cell,
-                dist.distance, std::abs(dist.surface_index));
+                dist.distance, std::abs(dist.surface));
 
               // Advance particle
               for (int lev = 0; lev < p.n_coord(); ++lev) {
                 p.coord(lev).r += dist.distance * p.coord(lev).u;
               }
-              p.surface() = dist.surface_index;
+              p.surface() = dist.surface;
               p.n_coord_last() = p.n_coord();
               p.n_coord() = dist.coord_level;
               if (dist.lattice_translation[0] != 0 ||

src/tallies/filter_musurface.cpp

@@ -12,7 +12,7 @@ void MuSurfaceFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
   // Get surface normal (and make sure it is a unit vector)
-  const auto surf {model::surfaces[std::abs(p.surface()) - 1].get()};
+  const auto surf {model::surfaces[p.surface_index()].get()};
   auto n = surf->normal(p.r());
   n /= n.norm();
 

src/tallies/filter_surface.cpp

@@ -47,7 +47,7 @@ void SurfaceFilter::set_surfaces(gsl::span<int32_t> surfaces)
 void SurfaceFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
-  auto search = map_.find(std::abs(p.surface()) - 1);
+  auto search = map_.find(p.surface_index());
   if (search != map_.end()) {
     match.bins_.push_back(search->second);
     if (p.surface() < 0) {