Tweaks to Universe::find_cell (#2662)

src/dagmc.cpp

@@ -182,10 +182,11 @@ void DAGUniverse::init_geometry()
       model::cell_map[c->id_] = model::cells.size();
     } else {
       warning(fmt::format("DAGMC Cell IDs: {}", dagmc_ids_for_dim(3)));
-      fatal_error(fmt::format("DAGMC Universe {} contains a cell with ID {}, which "
-                              "already exists elsewhere in the geometry. Setting auto_geom_ids "
-                              "to True when initiating the DAGMC Universe may "
-                              "resolve this issue",
+      fatal_error(fmt::format(
+        "DAGMC Universe {} contains a cell with ID {}, which "
+        "already exists elsewhere in the geometry. Setting auto_geom_ids "
+        "to True when initiating the DAGMC Universe may "
+        "resolve this issue",
         this->id_, c->id_));
     }
 
@@ -395,7 +396,7 @@ bool DAGUniverse::find_cell(Particle& p) const
   // cells, place it in the implicit complement
   bool found = Universe::find_cell(p);
   if (!found && model::universe_map[this->id_] != model::root_universe) {
-    p.coord(p.n_coord() - 1).cell = implicit_complement_idx();
+    p.lowest_coord().cell = implicit_complement_idx();
     found = true;
   }
   return found;
@@ -581,7 +582,7 @@ std::pair<double, int32_t> DAGCell::distance(
     p->history().reset();
   }
 
-  const auto& univ = model::universes[p->coord(p->n_coord() - 1).universe];
+  const auto& univ = model::universes[p->lowest_coord().universe];
 
   DAGUniverse* dag_univ = static_cast<DAGUniverse*>(univ.get());
   if (!dag_univ)
@@ -609,7 +610,8 @@ std::pair<double, int32_t> DAGCell::distance(
     // into the implicit complement on the other side where no intersection will
     // be found. Treating this as a lost particle is problematic when plotting.
     // Instead, the infinite distance and invalid surface index are returned.
-    if (settings::run_mode == RunMode::PLOTTING) return {INFTY, -1};
+    if (settings::run_mode == RunMode::PLOTTING)
+      return {INFTY, -1};
 
     // the particle should be marked as lost immediately if an intersection
     // isn't found in a volume that is not the implicit complement. In the case
@@ -739,7 +741,8 @@ void check_dagmc_root_univ()
   }
 }
 
-int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t 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 cellp = dynamic_cast<DAGCell*>(model::cells[curr_cell].get());
   auto univp = static_cast<DAGUniverse*>(model::universes[univ].get());
@@ -750,11 +753,12 @@ int32_t next_cell(int32_t surf, int32_t curr_cell, int32_t univ) {
     cellp->dagmc_ptr()->entity_by_index(3, cellp->dag_index());
 
   moab::EntityHandle new_vol;
-  moab::ErrorCode rval = cellp->dagmc_ptr()->next_vol(surf_handle, curr_vol, new_vol);
-  if (rval != moab::MB_SUCCESS) return -1;
+  moab::ErrorCode rval =
+    cellp->dagmc_ptr()->next_vol(surf_handle, curr_vol, new_vol);
+  if (rval != moab::MB_SUCCESS)
+    return -1;
 
-  return cellp->dagmc_ptr()->index_by_handle(new_vol) +
-         univp->cell_idx_offset_;
+  return cellp->dagmc_ptr()->index_by_handle(new_vol) + univp->cell_idx_offset_;
 }
 
 } // namespace openmc

src/geometry.cpp

@@ -110,7 +110,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list)
       i_cell = *it;
 
       // Make sure the search cell is in the same universe.
-      int i_universe = p.coord(p.n_coord() - 1).universe;
+      int i_universe = p.lowest_coord().universe;
       if (model::cells[i_cell]->universe_ != i_universe)
         continue;
 
@@ -119,7 +119,7 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list)
       Direction u {p.u_local()};
       auto surf = p.surface();
       if (model::cells[i_cell]->contains(r, u, surf)) {
-        p.coord(p.n_coord() - 1).cell = i_cell;
+        p.lowest_coord().cell = i_cell;
         found = true;
         break;
       }
@@ -145,15 +145,15 @@ bool find_cell_inner(Particle& p, const NeighborList* neighbor_list)
     // code below this conditional, we set i_cell back to C_NONE to indicate
     // that.
     if (i_cell == C_NONE) {
-      int i_universe = p.coord(p.n_coord() - 1).universe;
+      int i_universe = p.lowest_coord().universe;
       const auto& univ {model::universes[i_universe]};
       found = univ->find_cell(p);
     }
 
     if (!found) {
       return found;
     }
-    i_cell = p.coord(p.n_coord() - 1).cell;
+    i_cell = p.lowest_coord().cell;
 
     // Announce the cell that the particle is entering.
     if (found && (settings::verbosity >= 10 || p.trace())) {
@@ -289,7 +289,7 @@ bool neighbor_list_find_cell(Particle& p)
 
 bool exhaustive_find_cell(Particle& p)
 {
-  int i_universe = p.coord(p.n_coord() - 1).universe;
+  int i_universe = p.lowest_coord().universe;
   if (i_universe == C_NONE) {
     p.coord(0).universe = model::root_universe;
     p.n_coord() = 1;
@@ -306,7 +306,7 @@ bool exhaustive_find_cell(Particle& p)
 
 void cross_lattice(Particle& p, const BoundaryInfo& boundary)
 {
-  auto& coord {p.coord(p.n_coord() - 1)};
+  auto& coord {p.lowest_coord()};
   auto& lat {*model::lattices[coord.lattice]};
 
   if (settings::verbosity >= 10 || p.trace()) {
@@ -344,7 +344,7 @@ void cross_lattice(Particle& p, const BoundaryInfo& boundary)
 
   } else {
     // Find cell in next lattice element.
-    p.coord(p.n_coord() - 1).universe = lat[coord.lattice_i];
+    p.lowest_coord().universe = lat[coord.lattice_i];
     bool found = exhaustive_find_cell(p);
 
     if (!found) {
@@ -472,7 +472,7 @@ extern "C" int openmc_find_cell(
     return OPENMC_E_GEOMETRY;
   }
 
-  *index = p.coord(p.n_coord() - 1).cell;
+  *index = p.lowest_coord().cell;
   *instance = p.cell_instance();
   return 0;
 }

src/particle.cpp

@@ -147,7 +147,7 @@ void Particle::event_calculate_xs()
   // If the cell hasn't been determined based on the particle's location,
   // initiate a search for the current cell. This generally happens at the
   // beginning of the history and again for any secondary particles
-  if (coord(n_coord() - 1).cell == C_NONE) {
+  if (lowest_coord().cell == C_NONE) {
     if (!exhaustive_find_cell(*this)) {
       mark_as_lost(
         "Could not find the cell containing particle " + std::to_string(id()));
@@ -156,7 +156,7 @@ void Particle::event_calculate_xs()
 
     // Set birth cell attribute
     if (cell_born() == C_NONE)
-      cell_born() = coord(n_coord() - 1).cell;
+      cell_born() = lowest_coord().cell;
   }
 
   // Write particle track.
@@ -392,15 +392,15 @@ void Particle::event_revive_from_secondary()
       // Since the birth cell of the particle has not been set we
       // have to determine it before the energy of the secondary particle can be
       // removed from the pulse-height of this cell.
-      if (coord(n_coord() - 1).cell == C_NONE) {
+      if (lowest_coord().cell == C_NONE) {
         if (!exhaustive_find_cell(*this)) {
           mark_as_lost("Could not find the cell containing particle " +
                        std::to_string(id()));
           return;
         }
         // Set birth cell attribute
         if (cell_born() == C_NONE)
-          cell_born() = coord(n_coord() - 1).cell;
+          cell_born() = lowest_coord().cell;
       }
       pht_secondary_particles();
     }
@@ -456,7 +456,7 @@ void Particle::pht_collision_energy()
 
   // determine index of cell in pulse_height_cells
   auto it = std::find(model::pulse_height_cells.begin(),
-    model::pulse_height_cells.end(), coord(n_coord() - 1).cell);
+    model::pulse_height_cells.end(), lowest_coord().cell);
 
   if (it != model::pulse_height_cells.end()) {
     int index = std::distance(model::pulse_height_cells.begin(), it);
@@ -535,7 +535,7 @@ void Particle::cross_surface()
     material_last() = material();
     sqrtkT_last() = sqrtkT();
     // set new cell value
-    coord(n_coord() - 1).cell = i_cell;
+    lowest_coord().cell = i_cell;
     cell_instance() = 0;
     material() = model::cells[i_cell]->material_[0];
     sqrtkT() = model::cells[i_cell]->sqrtkT_[0];

src/plot.cpp

@@ -58,7 +58,7 @@ void IdData::set_value(size_t y, size_t x, const Particle& p, int level)
   }
 
   // set material data
-  Cell* c = model::cells.at(p.coord(p.n_coord() - 1).cell).get();
+  Cell* c = model::cells.at(p.lowest_coord().cell).get();
   if (p.material() == MATERIAL_VOID) {
     data_(y, x, 2) = MATERIAL_VOID;
     return;
@@ -79,7 +79,7 @@ PropertyData::PropertyData(size_t h_res, size_t v_res)
 
 void PropertyData::set_value(size_t y, size_t x, const Particle& p, int level)
 {
-  Cell* c = model::cells.at(p.coord(p.n_coord() - 1).cell).get();
+  Cell* c = model::cells.at(p.lowest_coord().cell).get();
   data_(y, x, 0) = (p.sqrtkT() * p.sqrtkT()) / K_BOLTZMANN;
   if (c->type_ != Fill::UNIVERSE && p.material() != MATERIAL_VOID) {
     Material* m = model::materials.at(p.material()).get();
@@ -1327,9 +1327,8 @@ void ProjectionPlot::create_output() const
               // edges on the model boundary for the same cell.
               if (first_inside_model) {
                 this_line_segments[tid][horiz].emplace_back(
-                  color_by_ == PlotColorBy::mats
-                    ? p.material()
-                    : p.coord(p.n_coord() - 1).cell,
+                  color_by_ == PlotColorBy::mats ? p.material()
+                                                 : p.lowest_coord().cell,
                   0.0, first_surface);
                 first_inside_model = false;
               }
@@ -1339,7 +1338,7 @@ void ProjectionPlot::create_output() const
               auto dist = distance_to_boundary(p);
               this_line_segments[tid][horiz].emplace_back(
                 color_by_ == PlotColorBy::mats ? p.material()
-                                               : p.coord(p.n_coord() - 1).cell,
+                                               : p.lowest_coord().cell,
                 dist.distance, std::abs(dist.surface_index));
 
               // Advance particle

src/tallies/filter_cell_instance.cpp

@@ -72,7 +72,7 @@ void CellInstanceFilter::set_cell_instances(gsl::span<CellInstance> instances)
 void CellInstanceFilter::get_all_bins(
   const Particle& p, TallyEstimator estimator, FilterMatch& match) const
 {
-  gsl::index index_cell = p.coord(p.n_coord() - 1).cell;
+  gsl::index index_cell = p.lowest_coord().cell;
   gsl::index instance = p.cell_instance();
 
   if (cells_.count(index_cell) > 0) {

src/universe.cpp

@@ -41,18 +41,17 @@ bool Universe::find_cell(Particle& p) const
   const auto& cells {
     !partitioner_ ? cells_ : partitioner_->get_cells(p.r_local(), p.u_local())};
 
-  for (auto it = cells.begin(); it != cells.end(); it++) {
-    int32_t i_cell = *it;
-    int32_t i_univ = p.coord(p.n_coord() - 1).universe;
+  Position r {p.r_local()};
+  Position u {p.u_local()};
+  auto surf = p.surface();
+  int32_t i_univ = p.lowest_coord().universe;
+
+  for (auto i_cell : cells) {
     if (model::cells[i_cell]->universe_ != i_univ)
       continue;
-
-    // Check if this cell contains the particle;
-    Position r {p.r_local()};
-    Direction u {p.u_local()};
-    auto surf = p.surface();
+    // Check if this cell contains the particle
     if (model::cells[i_cell]->contains(r, u, surf)) {
-      p.coord(p.n_coord() - 1).cell = i_cell;
+      p.lowest_coord().cell = i_cell;
       return true;
     }
   }