From d0a0d34b89cc59329d60bf26b657768a7ccd64e3 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Wed, 3 Apr 2019 19:50:12 -0400 Subject: [PATCH] Address #1210 comments --- include/openmc/cell.h | 37 +---------------- src/cell.cpp | 96 ++++++++++++++++++++++++++++++++----------- src/geometry_aux.cpp | 5 +-- 3 files changed, 76 insertions(+), 62 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 2cd0bd3fa0..cc4429b44f 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -210,44 +210,9 @@ public: explicit UniversePartitioner(const Universe& univ); //! Return the list of cells that could contain the given coordinates. - const std::vector& get_cells(Position r, Direction u) const - {return get_cells(r, u, 0, surfs_.size()-1, (surfs_.size()-1)/2);} + const std::vector& get_cells(Position r, Direction u) const; private: - //! Perform a binary search for the given coordinates. - // - //! \param left The lower search bound of the `surfs_` vector - //! \param right The upper search bound of the `surfs_` vector - //! \param middle The index of the `Surface` to test against - const std::vector& - get_cells(Position r, Direction u, int left, int right, int middle) const - { - // Check the sense of the coordinates for the current surface. - auto i_surf = surfs_[middle]; - const auto& surf = *model::surfaces[i_surf]; - if (surf.sense(r, u)) { - // The coordinates lie in the positive halfspace. Recurse if there are - // more surfaces to check. Otherwise, return the cells on the positive - // side of this surface. - int right_leaf = right - (right - middle) / 2; - if (right_leaf != middle) { - return get_cells(r, u, middle+1, right, right_leaf); - } else { - return partitions_[middle+1]; - } - } else { - // The coordinates lie in the negative halfspace. Recurse if there are - // more surfaces to check. Otherwise, return the cells on the negative - // side of this surface. - int left_leaf = left + (middle - left) / 2; - if (left_leaf != middle) { - return get_cells(r, u, left, middle-1, left_leaf); - } else { - return partitions_[middle]; - } - } - } - //! A sorted vector of indices to surfaces that partition the universe std::vector surfs_; diff --git a/src/cell.cpp b/src/cell.cpp index 5e803c4cd5..8e2dfc5a7a 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -659,10 +659,10 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) } } - // Define a functor for comparing z-planes by their location, and use it to + // Define a lambda for comparing z-planes by their location, and use it to // sort the surfs_ member. - struct { - bool operator()(int32_t i_surf, int32_t j_surf) const + std::sort(surfs_.begin(), surfs_.end(), + [](int32_t i_surf, int32_t j_surf) -> bool { const auto* surf = model::surfaces[i_surf].get(); const auto* zplane = dynamic_cast(surf); @@ -672,25 +672,24 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) double zj = zplane->z0_; return zi < zj; } - } compare_surfs; - std::sort(surfs_.begin(), surfs_.end(), compare_surfs); + ); // Populate the partition lists. partitions_.resize(surfs_.size() + 1); for (auto i_cell : univ.cells_) { // Find the tokens for bounding z-planes. - int32_t min_token = 0, max_token = 0; + int32_t lower_token = 0, upper_token = 0; double min_z, max_z; for (auto token : model::cells[i_cell]->rpn_) { if (token < OP_UNION) { const auto* surf = model::surfaces[std::abs(token) - 1].get(); if (const auto* zplane = dynamic_cast(surf)) { - if (min_token == 0 || zplane->z0_ < min_z) { - min_token = token; + if (lower_token == 0 || zplane->z0_ < min_z) { + lower_token = token; min_z = zplane->z0_; } - if (max_token == 0 || zplane->z0_ > max_z) { - max_token = token; + if (upper_token == 0 || zplane->z0_ > max_z) { + upper_token = token; max_z = zplane->z0_; } } @@ -698,27 +697,78 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) } // If there are no bounding z-planes, add this cell to all partitions. - if (min_token == 0) { + if (lower_token == 0) { for (auto& p : partitions_) p.push_back(i_cell); continue; } - // Iterate over partitions, and add this cell to each appropriate one. - // Since surfs_ is sorted, we know this cell belongs to every partition - // between min_token and max_token. - bool in_partition = min_token < 0; - for (auto i = 0; i < surfs_.size(); ++i) { - if (in_partition) { - partitions_[i].push_back(i_cell); - if (max_token == -(surfs_[i] + 1)) { - in_partition = false; + // Find the first partition this cell lies in. If the lower_token indicates + // a negative halfspace, then the cell is unbounded in the lower direction + // and it lies in the first partition onward. Otherwise, it is bounded by + // the positive halfspace given by the lower_token. + int first_partition = 0; + if (lower_token > 0) { + for (int i = 0; i < surfs_.size(); ++i) { + if (lower_token == surfs_[i] + 1) { + first_partition = i + 1; break; } - } else { - in_partition = min_token == surfs_[i] + 1; } } - if (in_partition) partitions_.back().push_back(i_cell); + + // Find the last partition this cell lies in. The logic is analogous to the + // logic for first_partition. + int last_partition = surfs_.size(); + if (upper_token < 0) { + for (int i = first_partition; i < surfs_.size(); ++i) { + if (upper_token == -(surfs_[i] + 1)) { + last_partition = i; + break; + } + } + } + + // Add the cell to all relevant partitions. + for (int i = first_partition; i < last_partition + 1; ++i) { + partitions_[i].push_back(i_cell); + } + } +} + +const std::vector& +UniversePartitioner::get_cells(Position r, Direction u) const +{ + // Perform a binary search for the partition containing the given coordinates. + int left = 0; + int middle = (surfs_.size() - 1) / 2; + int right = surfs_.size() - 1; + while (true) { + // Check the sense of the coordinates for the current surface. + const auto& surf = *model::surfaces[surfs_[middle]]; + if (surf.sense(r, u)) { + // The coordinates lie in the positive halfspace. Recurse if there are + // more surfaces to check. Otherwise, return the cells on the positive + // side of this surface. + int right_leaf = right - (right - middle) / 2; + if (right_leaf != middle) { + left = middle + 1; + middle = right_leaf; + } else { + return partitions_[middle+1]; + } + + } else { + // The coordinates lie in the negative halfspace. Recurse if there are + // more surfaces to check. Otherwise, return the cells on the negative + // side of this surface. + int left_leaf = left + (middle - left) / 2; + if (left_leaf != middle) { + right = middle-1; + middle = left_leaf; + } else { + return partitions_[middle]; + } + } } } diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 38d3a8a0e7..f02ee42b14 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -139,14 +139,13 @@ partition_universes() } // Partition the universe if there are more than 5 z-planes. (Fewer than - // five is likely no worth it.) + // 5 is likely not worth it.) int n_zplanes = 0; for (auto i_surf : surf_inds) { if (dynamic_cast(model::surfaces[i_surf].get())) { ++n_zplanes; if (n_zplanes > 5) { - univ->partitioner_ = std::make_unique( - UniversePartitioner(*univ)); + univ->partitioner_ = std::make_unique(*univ); break; } }