diff --git a/CMakeLists.txt b/CMakeLists.txt index e972c49378..788e8c6814 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -351,6 +351,7 @@ list(APPEND libopenmc_SOURCES src/timer.cpp src/thermal.cpp src/track_output.cpp + src/universe.cpp src/urr.cpp src/volume_calc.cpp src/weight_windows.cpp diff --git a/include/openmc/cell.h b/include/openmc/cell.h index afc9a33c59..d975750a16 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -17,6 +17,7 @@ #include "openmc/neighbor_list.h" #include "openmc/position.h" #include "openmc/surface.h" +#include "openmc/universe.h" #include "openmc/vector.h" namespace openmc { @@ -48,36 +49,8 @@ namespace model { extern std::unordered_map cell_map; extern vector> cells; -extern std::unordered_map universe_map; -extern vector> universes; } // namespace model -//============================================================================== -//! A geometry primitive that fills all space and contains cells. -//============================================================================== - -class Universe { -public: - int32_t id_; //!< Unique ID - vector cells_; //!< Cells within this universe - - //! \brief Write universe information to an HDF5 group. - //! \param group_id An HDF5 group id. - virtual void to_hdf5(hid_t group_id) const; - - virtual bool find_cell(Particle& p) const; - - BoundingBox bounding_box() const; - - const GeometryType& geom_type() const { return geom_type_; } - GeometryType& geom_type() { return geom_type_; } - - unique_ptr partitioner_; - -private: - GeometryType geom_type_ = GeometryType::CSG; -}; - //============================================================================== //============================================================================== @@ -296,35 +269,6 @@ protected: vector::iterator start, const vector& rpn); }; -//============================================================================== -//! Speeds up geometry searches by grouping cells in a search tree. -// -//! Currently this object only works with universes that are divided up by a -//! bunch of z-planes. It could be generalized to other planes, cylinders, -//! and spheres. -//============================================================================== - -class UniversePartitioner { -public: - explicit UniversePartitioner(const Universe& univ); - - //! Return the list of cells that could contain the given coordinates. - const vector& get_cells(Position r, Direction u) const; - -private: - //! A sorted vector of indices to surfaces that partition the universe - vector surfs_; - - //! Vectors listing the indices of the cells that lie within each partition - // - //! There are n+1 partitions with n surfaces. `partitions_.front()` gives the - //! cells that lie on the negative side of `surfs_.front()`. - //! `partitions_.back()` gives the cells that lie on the positive side of - //! `surfs_.back()`. Otherwise, `partitions_[i]` gives cells sandwiched - //! between `surfs_[i-1]` and `surfs_[i]`. - vector> partitions_; -}; - //============================================================================== //! Define an instance of a particular cell //============================================================================== @@ -356,9 +300,5 @@ struct CellInstanceHash { void read_cells(pugi::xml_node node); -#ifdef DAGMC -class DAGUniverse; -#endif - } // namespace openmc #endif // OPENMC_CELL_H diff --git a/include/openmc/universe.h b/include/openmc/universe.h new file mode 100644 index 0000000000..b6ed3e7246 --- /dev/null +++ b/include/openmc/universe.h @@ -0,0 +1,79 @@ +#ifndef OPENMC_UNIVERSE_H +#define OPENMC_UNIVERSE_H + +#include "openmc/cell.h" + +namespace openmc { + +#ifdef DAGMC +class DAGUniverse; +#endif + +class Universe; +class UniversePartitioner; + +namespace model { + +extern std::unordered_map universe_map; +extern vector> universes; + +} // namespace model + + +//============================================================================== +//! A geometry primitive that fills all space and contains cells. +//============================================================================== + +class Universe { +public: + int32_t id_; //!< Unique ID + vector cells_; //!< Cells within this universe + + //! \brief Write universe information to an HDF5 group. + //! \param group_id An HDF5 group id. + virtual void to_hdf5(hid_t group_id) const; + + virtual bool find_cell(Particle& p) const; + + BoundingBox bounding_box() const; + + const GeometryType& geom_type() const { return geom_type_; } + GeometryType& geom_type() { return geom_type_; } + + unique_ptr partitioner_; + +private: + GeometryType geom_type_ = GeometryType::CSG; +}; + +//============================================================================== +//! Speeds up geometry searches by grouping cells in a search tree. +// +//! Currently this object only works with universes that are divided up by a +//! bunch of z-planes. It could be generalized to other planes, cylinders, +//! and spheres. +//============================================================================== + +class UniversePartitioner { +public: + explicit UniversePartitioner(const Universe& univ); + + //! Return the list of cells that could contain the given coordinates. + const vector& get_cells(Position r, Direction u) const; + +private: + //! A sorted vector of indices to surfaces that partition the universe + vector surfs_; + + //! Vectors listing the indices of the cells that lie within each partition + // + //! There are n+1 partitions with n surfaces. `partitions_.front()` gives the + //! cells that lie on the negative side of `surfs_.front()`. + //! `partitions_.back()` gives the cells that lie on the positive side of + //! `surfs_.back()`. Otherwise, `partitions_[i]` gives cells sandwiched + //! between `surfs_[i-1]` and `surfs_[i]`. + vector> partitions_; +}; + +} // namespace openmc +#endif // OPENMC_UNIVERSE_H \ No newline at end of file diff --git a/src/cell.cpp b/src/cell.cpp index 2b6555d128..28be2c2cff 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -34,8 +34,6 @@ namespace model { std::unordered_map cell_map; vector> cells; -std::unordered_map universe_map; -vector> universes; } // namespace model //============================================================================== @@ -186,66 +184,6 @@ vector generate_rpn(int32_t cell_id, vector infix) return rpn; } -//============================================================================== -// Universe implementation -//============================================================================== - -void Universe::to_hdf5(hid_t universes_group) const -{ - // Create a group for this universe. - auto group = create_group(universes_group, fmt::format("universe {}", id_)); - - // Write the geometry representation type. - write_string(group, "geom_type", "csg", false); - - // Write the contained cells. - if (cells_.size() > 0) { - vector cell_ids; - for (auto i_cell : cells_) - cell_ids.push_back(model::cells[i_cell]->id_); - write_dataset(group, "cells", cell_ids); - } - - close_group(group); -} - -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; - 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(); - if (model::cells[i_cell]->contains(r, u, surf)) { - p.coord(p.n_coord() - 1).cell = i_cell; - return true; - } - } - return false; -} - -BoundingBox Universe::bounding_box() const -{ - BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY}; - if (cells_.size() == 0) { - return {}; - } else { - for (const auto& cell : cells_) { - auto& c = model::cells[cell]; - bbox |= c->bounding_box(); - } - } - return bbox; -} - //============================================================================== // Cell implementation //============================================================================== @@ -876,151 +814,6 @@ bool CSGCell::contains_complex( } } -//============================================================================== -// UniversePartitioner implementation -//============================================================================== - -UniversePartitioner::UniversePartitioner(const Universe& univ) -{ - // Define an ordered set of surface indices that point to z-planes. Use a - // functor to to order the set by the z0_ values of the corresponding planes. - struct compare_surfs { - bool operator()(const int32_t& i_surf, const int32_t& j_surf) const - { - const auto* surf = model::surfaces[i_surf].get(); - const auto* zplane = dynamic_cast(surf); - double zi = zplane->z0_; - surf = model::surfaces[j_surf].get(); - zplane = dynamic_cast(surf); - double zj = zplane->z0_; - return zi < zj; - } - }; - std::set surf_set; - - // Find all of the z-planes in this universe. A set is used here for the - // O(log(n)) insertions that will ensure entries are not repeated. - for (auto i_cell : univ.cells_) { - for (auto token : model::cells[i_cell]->rpn_) { - if (token < OP_UNION) { - auto i_surf = std::abs(token) - 1; - const auto* surf = model::surfaces[i_surf].get(); - if (const auto* zplane = dynamic_cast(surf)) - surf_set.insert(i_surf); - } - } - } - - // Populate the surfs_ vector from the ordered set. - surfs_.insert(surfs_.begin(), surf_set.begin(), surf_set.end()); - - // Populate the partition lists. - partitions_.resize(surfs_.size() + 1); - for (auto i_cell : univ.cells_) { - // It is difficult to determine the bounds of a complex cell, so add complex - // cells to all partitions. - if (!model::cells[i_cell]->simple_) { - for (auto& p : partitions_) - p.push_back(i_cell); - continue; - } - - // Find the tokens for bounding z-planes. - 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 (lower_token == 0 || zplane->z0_ < min_z) { - lower_token = token; - min_z = zplane->z0_; - } - if (upper_token == 0 || zplane->z0_ > max_z) { - upper_token = token; - max_z = zplane->z0_; - } - } - } - } - - // If there are no bounding z-planes, add this cell to all partitions. - if (lower_token == 0) { - for (auto& p : partitions_) - p.push_back(i_cell); - continue; - } - - // 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; - } - } - } - - // 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; ++i) { - partitions_[i].push_back(i_cell); - } - } -} - -const 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]; - } - } - } -} - //============================================================================== // Non-method functions //============================================================================== diff --git a/src/universe.cpp b/src/universe.cpp new file mode 100644 index 0000000000..f4b5e7d29b --- /dev/null +++ b/src/universe.cpp @@ -0,0 +1,221 @@ +#include "openmc/universe.h" + +#include + +#include "openmc/hdf5_interface.h" + +namespace openmc { + +namespace model { + +std::unordered_map universe_map; +vector> universes; + +} // namespace model + +//============================================================================== +// Universe implementation +//============================================================================== + +void Universe::to_hdf5(hid_t universes_group) const +{ + // Create a group for this universe. + auto group = create_group(universes_group, fmt::format("universe {}", id_)); + + // Write the geometry representation type. + write_string(group, "geom_type", "csg", false); + + // Write the contained cells. + if (cells_.size() > 0) { + vector cell_ids; + for (auto i_cell : cells_) + cell_ids.push_back(model::cells[i_cell]->id_); + write_dataset(group, "cells", cell_ids); + } + + close_group(group); +} + +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; + 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(); + if (model::cells[i_cell]->contains(r, u, surf)) { + p.coord(p.n_coord() - 1).cell = i_cell; + return true; + } + } + return false; +} + +BoundingBox Universe::bounding_box() const +{ + BoundingBox bbox = {INFTY, -INFTY, INFTY, -INFTY, INFTY, -INFTY}; + if (cells_.size() == 0) { + return {}; + } else { + for (const auto& cell : cells_) { + auto& c = model::cells[cell]; + bbox |= c->bounding_box(); + } + } + return bbox; +} + +//============================================================================== +// UniversePartitioner implementation +//============================================================================== + +UniversePartitioner::UniversePartitioner(const Universe& univ) +{ + // Define an ordered set of surface indices that point to z-planes. Use a + // functor to to order the set by the z0_ values of the corresponding planes. + struct compare_surfs { + bool operator()(const int32_t& i_surf, const int32_t& j_surf) const + { + const auto* surf = model::surfaces[i_surf].get(); + const auto* zplane = dynamic_cast(surf); + double zi = zplane->z0_; + surf = model::surfaces[j_surf].get(); + zplane = dynamic_cast(surf); + double zj = zplane->z0_; + return zi < zj; + } + }; + std::set surf_set; + + // Find all of the z-planes in this universe. A set is used here for the + // O(log(n)) insertions that will ensure entries are not repeated. + for (auto i_cell : univ.cells_) { + for (auto token : model::cells[i_cell]->rpn_) { + if (token < OP_UNION) { + auto i_surf = std::abs(token) - 1; + const auto* surf = model::surfaces[i_surf].get(); + if (const auto* zplane = dynamic_cast(surf)) + surf_set.insert(i_surf); + } + } + } + + // Populate the surfs_ vector from the ordered set. + surfs_.insert(surfs_.begin(), surf_set.begin(), surf_set.end()); + + // Populate the partition lists. + partitions_.resize(surfs_.size() + 1); + for (auto i_cell : univ.cells_) { + // It is difficult to determine the bounds of a complex cell, so add complex + // cells to all partitions. + if (!model::cells[i_cell]->simple_) { + for (auto& p : partitions_) + p.push_back(i_cell); + continue; + } + + // Find the tokens for bounding z-planes. + 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 (lower_token == 0 || zplane->z0_ < min_z) { + lower_token = token; + min_z = zplane->z0_; + } + if (upper_token == 0 || zplane->z0_ > max_z) { + upper_token = token; + max_z = zplane->z0_; + } + } + } + } + + // If there are no bounding z-planes, add this cell to all partitions. + if (lower_token == 0) { + for (auto& p : partitions_) + p.push_back(i_cell); + continue; + } + + // 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; + } + } + } + + // 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; ++i) { + partitions_[i].push_back(i_cell); + } + } +} + +const 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]; + } + } + } +} + +} \ No newline at end of file