From 9127ff4f5ecacd46692784f1a75e490490d1b3f2 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Thu, 8 Sep 2022 11:24:08 -0500 Subject: [PATCH 1/5] [WIP] started region refactor --- include/openmc/cell.h | 40 ++++++++++++++++++++++++++++++++++++++-- src/cell.cpp | 7 +++++-- 2 files changed, 43 insertions(+), 4 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 8171967273..fcce520fed 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -52,8 +52,41 @@ extern vector> cells; } // namespace model //============================================================================== + +// TODO: Maybe not, just move this inline to Region::bounding_box() for complex +// case +class RegionPostfix { +public: + BoundingBox bounding_box() const; +}; + +class Region { +public: + Region() {} + explicit Region(std::string region_expression); + + void add_precedence(); + std::string str() const; + BoundingBox bounding_box() const; + bool contains(Position r, Direction u, int32_t on_surface) const; + + RegionPostfix to_postfix() const; + +private: + void add_parentheses(); + void apply_demorgan( + vector::iterator start, vector::iterator stop); + + //! Definition of spatial region as Boolean expression of half-spaces + // TODO: Should this be a vector of some other type + vector tokens_; + bool simple_; //!< Does the region contain only intersections? +}; + //============================================================================== +// TODO: Think about what data members really need to live in this class versus +// putting them in the Region class class Cell { public: //---------------------------------------------------------------------------- @@ -198,9 +231,9 @@ public: //! T. The units are sqrt(eV). vector sqrtkT_; + // TODO: Probably move this guy to CSGCell //! Definition of spatial region as Boolean expression of half-spaces - vector region_; - bool simple_; //!< Does the region contain only intersections? + Region region_; //! \brief Neighboring cells in the same universe. NeighborList neighbors_; @@ -263,6 +296,9 @@ protected: //! \param rpn The rpn being searched static vector::iterator find_left_parenthesis( vector::iterator start, const vector& rpn); + +private: + Region region_; }; //============================================================================== diff --git a/src/cell.cpp b/src/cell.cpp index 4ff35498c6..4eb0fa0b95 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -43,6 +43,7 @@ vector> cells; //! operators. //============================================================================== +// TODO: Move this to be Region::Region(...) vector tokenize(const std::string region_spec) { // Check for an empty region_spec first. @@ -183,7 +184,8 @@ void add_precedence(std::vector& infix) // Set the current operator if is hasn't been set current_op = token; } else if (token != current_op) { - // If the current operator doesn't match the token, add parenthesis to assert precedence + // If the current operator doesn't match the token, add parenthesis to + // assert precedence it = add_parentheses(it, infix); current_op = 0; } @@ -586,7 +588,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // Get a tokenized representation of the region specification and apply De // Morgans law - region_ = tokenize(region_spec); + region_ = Region(region_spec); remove_complement_ops(region_); // Convert user IDs to surface indices. @@ -696,6 +698,7 @@ void CSGCell::to_hdf5_inner(hid_t group_id) const write_string(group_id, "geom_type", "csg", false); + // TODO: Move to Region::str() // Write the region specification. if (!region_.empty()) { std::stringstream region_spec {}; From 06312b0535bbc0d46464a77039aa0ad483c9db21 Mon Sep 17 00:00:00 2001 From: myerspat Date: Sun, 25 Sep 2022 22:05:22 -0400 Subject: [PATCH 2/5] Refactored CSGCell to use Region --- include/openmc/cell.h | 90 +++-- src/cell.cpp | 872 ++++++++++++++++++++++-------------------- src/geometry.cpp | 2 +- src/geometry_aux.cpp | 5 +- src/universe.cpp | 36 +- 5 files changed, 523 insertions(+), 482 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index fcce520fed..0891c062f3 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -53,33 +53,46 @@ extern vector> cells; //============================================================================== -// TODO: Maybe not, just move this inline to Region::bounding_box() for complex -// case -class RegionPostfix { -public: - BoundingBox bounding_box() const; -}; - class Region { public: + //---------------------------------------------------------------------------- + // Constructors Region() {} - explicit Region(std::string region_expression); + explicit Region(std::string region_expressioni, int32_t cell_id); - void add_precedence(); - std::string str() const; - BoundingBox bounding_box() const; + //---------------------------------------------------------------------------- + // Methods + std::pair distance( + Position r, Direction u, int32_t on_surface, Particle* p) const; bool contains(Position r, Direction u, int32_t on_surface) const; + BoundingBox bounding_box(int32_t cell_id) const; + std::string str() const; + vector surfaces() const; - RegionPostfix to_postfix() const; + //---------------------------------------------------------------------------- + // Accessors + bool is_simple() const { return simple_; } private: - void add_parentheses(); + //---------------------------------------------------------------------------- + // Private Methods + vector generate_postfix(int32_t cell_id) const; + bool contains_simple(Position r, Direction u, int32_t on_surface) const; + bool contains_complex(Position r, Direction u, int32_t on_surface) const; + BoundingBox bounding_box_simple() const; + BoundingBox bounding_box_complex(vector postfix) const; + void add_precedence(); + std::vector::iterator add_parentheses( + std::vector::iterator start); + void remove_complement_ops(); void apply_demorgan( vector::iterator start, vector::iterator stop); + //---------------------------------------------------------------------------- + // Private Data //! Definition of spatial region as Boolean expression of half-spaces // TODO: Should this be a vector of some other type - vector tokens_; + vector expression_; bool simple_; //!< Does the region contain only intersections? }; @@ -141,6 +154,12 @@ public: //! Get the BoundingBox for this cell. virtual BoundingBox bounding_box() const = 0; + //! Get a vector of surfaces in the cell + virtual vector surfaces() const = 0; + + //! Check if the cell region expression is simple + virtual bool is_simple() const = 0; + //---------------------------------------------------------------------------- // Accessors @@ -231,10 +250,6 @@ public: //! T. The units are sqrt(eV). vector sqrtkT_; - // TODO: Probably move this guy to CSGCell - //! Definition of spatial region as Boolean expression of half-spaces - Region region_; - //! \brief Neighboring cells in the same universe. NeighborList neighbors_; @@ -260,35 +275,36 @@ struct CellInstanceItem { class CSGCell : public Cell { public: + //---------------------------------------------------------------------------- + // Constructors CSGCell(); - explicit CSGCell(pugi::xml_node cell_node); - bool contains(Position r, Direction u, int32_t on_surface) const override; + //---------------------------------------------------------------------------- + // Methods + vector surfaces() const override { return region_.surfaces(); } std::pair distance( - Position r, Direction u, int32_t on_surface, Particle* p) const override; + Position r, Direction u, int32_t on_surface, Particle* p) const override + { + return region_.distance(r, u, on_surface, p); + } + + bool contains(Position r, Direction u, int32_t on_surface) const override + { + return region_.contains(r, u, on_surface); + } + + BoundingBox bounding_box() const override + { + return region_.bounding_box(id_); + } void to_hdf5_inner(hid_t group_id) const override; - BoundingBox bounding_box() const override; + bool is_simple() const override { return region_.is_simple(); } protected: - bool contains_simple(Position r, Direction u, int32_t on_surface) const; - bool contains_complex(Position r, Direction u, int32_t on_surface) const; - BoundingBox bounding_box_simple() const; - static BoundingBox bounding_box_complex(vector postfix); - - //! Applies DeMorgan's laws to a section of the RPN - //! \param start Starting point for token modification - //! \param stop Stopping point for token modification - static void apply_demorgan( - vector::iterator start, vector::iterator stop); - - //! Removes complement operators from the RPN - //! \param rpn The rpn to remove complement operators from. - static void remove_complement_ops(vector& rpn); - //! Returns the beginning position of a parenthesis block (immediately before //! two surface tokens) in the RPN given a starting position at the end of //! that block (immediately after two surface tokens) diff --git a/src/cell.cpp b/src/cell.cpp index 4eb0fa0b95..dc38a567ee 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -36,243 +36,6 @@ vector> cells; } // namespace model -//============================================================================== -//! Convert region specification string to integer tokens. -//! -//! The characters (, ), |, and ~ count as separate tokens since they represent -//! operators. -//============================================================================== - -// TODO: Move this to be Region::Region(...) -vector tokenize(const std::string region_spec) -{ - // Check for an empty region_spec first. - vector tokens; - if (region_spec.empty()) { - return tokens; - } - - // Parse all halfspaces and operators except for intersection (whitespace). - for (int i = 0; i < region_spec.size();) { - if (region_spec[i] == '(') { - tokens.push_back(OP_LEFT_PAREN); - i++; - - } else if (region_spec[i] == ')') { - tokens.push_back(OP_RIGHT_PAREN); - i++; - - } else if (region_spec[i] == '|') { - tokens.push_back(OP_UNION); - i++; - - } else if (region_spec[i] == '~') { - tokens.push_back(OP_COMPLEMENT); - i++; - - } else if (region_spec[i] == '-' || region_spec[i] == '+' || - std::isdigit(region_spec[i])) { - // This is the start of a halfspace specification. Iterate j until we - // find the end, then push-back everything between i and j. - int j = i + 1; - while (j < region_spec.size() && std::isdigit(region_spec[j])) { - j++; - } - tokens.push_back(std::stoi(region_spec.substr(i, j - i))); - i = j; - - } else if (std::isspace(region_spec[i])) { - i++; - - } else { - auto err_msg = - fmt::format("Region specification contains invalid character, \"{}\"", - region_spec[i]); - fatal_error(err_msg); - } - } - - // Add in intersection operators where a missing operator is needed. - int i = 0; - while (i < tokens.size() - 1) { - bool left_compat {(tokens[i] < OP_UNION) || (tokens[i] == OP_RIGHT_PAREN)}; - bool right_compat {(tokens[i + 1] < OP_UNION) || - (tokens[i + 1] == OP_LEFT_PAREN) || - (tokens[i + 1] == OP_COMPLEMENT)}; - if (left_compat && right_compat) { - tokens.insert(tokens.begin() + i + 1, OP_INTERSECTION); - } - i++; - } - - return tokens; -} - -//============================================================================== -//! Add precedence for infix regions so intersections have higher -//! precedence than unions using parentheses. -//============================================================================== - -std::vector::iterator add_parentheses( - std::vector::iterator start, std::vector& infix) -{ - int32_t start_token = *start; - // Add left parenthesis - if (start_token == OP_INTERSECTION) { - start = infix.insert(start - 1, OP_LEFT_PAREN); - } else { - start = infix.insert(start + 1, OP_LEFT_PAREN); - } - start++; - - // Initialize return iterator - auto return_iterator = infix.begin(); - - // Add right parenthesis - // While the start iterator is within the bounds of infix - while (start < infix.end()) { - start++; - - // If the current token is an operator and is different than the start token - if (*start >= OP_UNION && *start != start_token) { - // Skip wrapped regions but save iterator position to check precedence and - // add right parenthesis, right parenthesis position depends on the - // operator, when the operator is a union then do not include the operator - // in the region, when the operator is an intersection then include the - // operato and next surface - if (*start == OP_LEFT_PAREN) { - return_iterator = start; - int depth = 1; - do { - start++; - if (*start > OP_COMPLEMENT) { - if (*start == OP_RIGHT_PAREN) { - depth--; - } else { - depth++; - } - } - } while (depth > 0); - } else { - start = infix.insert( - start_token == OP_UNION ? start - 1 : start, OP_RIGHT_PAREN); - if (return_iterator == infix.begin()) { - return_iterator = start - 1; - } - return return_iterator; - } - } - } - // If we get here a right parenthesis hasn't been placed, - // return iterator - infix.push_back(OP_RIGHT_PAREN); - if (return_iterator == infix.begin()) { - return_iterator = start - 1; - } - return return_iterator; -} - -void add_precedence(std::vector& infix) -{ - int32_t current_op = 0; - - for (auto it = infix.begin(); it != infix.end(); it++) { - int32_t token = *it; - - if (token == OP_UNION || token == OP_INTERSECTION) { - if (current_op == 0) { - // Set the current operator if is hasn't been set - current_op = token; - } else if (token != current_op) { - // If the current operator doesn't match the token, add parenthesis to - // assert precedence - it = add_parentheses(it, infix); - current_op = 0; - } - } else if (token > OP_COMPLEMENT) { - // If the token is a parenthesis reset the current operator - current_op = 0; - } - } -} - -//============================================================================== -//! Convert infix region specification to Reverse Polish Notation (RPN) -//! -//! This function uses the shunting-yard algorithm. -//============================================================================== - -vector generate_postfix(int32_t cell_id, vector infix) -{ - vector rpn; - vector stack; - - for (int32_t token : infix) { - if (token < OP_UNION) { - // If token is not an operator, add it to output - rpn.push_back(token); - } else if (token < OP_RIGHT_PAREN) { - // Regular operators union, intersection, complement - while (stack.size() > 0) { - int32_t op = stack.back(); - - if (op < OP_RIGHT_PAREN && ((token == OP_COMPLEMENT && token < op) || - (token != OP_COMPLEMENT && token <= op))) { - // While there is an operator, op, on top of the stack, if the token - // is left-associative and its precedence is less than or equal to - // that of op or if the token is right-associative and its precedence - // is less than that of op, move op to the output queue and push the - // token on to the stack. Note that only complement is - // right-associative. - rpn.push_back(op); - stack.pop_back(); - } else { - break; - } - } - - stack.push_back(token); - - } else if (token == OP_LEFT_PAREN) { - // If the token is a left parenthesis, push it onto the stack - stack.push_back(token); - - } else { - // If the token is a right parenthesis, move operators from the stack to - // the output queue until reaching the left parenthesis. - for (auto it = stack.rbegin(); *it != OP_LEFT_PAREN; it++) { - // If we run out of operators without finding a left parenthesis, it - // means there are mismatched parentheses. - if (it == stack.rend()) { - fatal_error(fmt::format( - "Mismatched parentheses in region specification for cell {}", - cell_id)); - } - rpn.push_back(stack.back()); - stack.pop_back(); - } - - // Pop the left parenthesis. - stack.pop_back(); - } - } - - while (stack.size() > 0) { - int32_t op = stack.back(); - - // If the operator is a parenthesis it is mismatched. - if (op >= OP_RIGHT_PAREN) { - fatal_error(fmt::format( - "Mismatched parentheses in region specification for cell {}", cell_id)); - } - - rpn.push_back(stack.back()); - stack.pop_back(); - } - - return rpn; -} - //============================================================================== // Cell implementation //============================================================================== @@ -588,44 +351,8 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // Get a tokenized representation of the region specification and apply De // Morgans law - region_ = Region(region_spec); - remove_complement_ops(region_); - - // Convert user IDs to surface indices. - for (auto& r : region_) { - if (r < OP_UNION) { - const auto& it {model::surface_map.find(abs(r))}; - if (it == model::surface_map.end()) { - throw std::runtime_error { - "Invalid surface ID " + std::to_string(abs(r)) + - " specified in region for cell " + std::to_string(id_) + "."}; - } - r = (r > 0) ? it->second + 1 : -(it->second + 1); - } - } - - // Check if this is a simple cell. - simple_ = true; - for (int32_t token : region_) { - if (token == OP_UNION) { - simple_ = false; - // Ensure intersections have precedence over unions - add_precedence(region_); - break; - } - } - region_.shrink_to_fit(); - - // If this cell is simple, remove all the superfluous operator tokens. - if (simple_) { - for (auto it = region_.begin(); it != region_.end(); it++) { - if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) { - region_.erase(it); - it--; - } - } - region_.shrink_to_fit(); - } + Region region(region_spec, id_); + region_ = region; // Read the translation vector. if (check_for_node(cell_node, "translation")) { @@ -652,99 +379,13 @@ CSGCell::CSGCell(pugi::xml_node cell_node) //============================================================================== -bool CSGCell::contains(Position r, Direction u, int32_t on_surface) const -{ - if (simple_) { - return contains_simple(r, u, on_surface); - } else { - return contains_complex(r, u, on_surface); - } -} - -//============================================================================== - -std::pair CSGCell::distance( - Position r, Direction u, int32_t on_surface, Particle* p) const -{ - double min_dist {INFTY}; - int32_t i_surf {std::numeric_limits::max()}; - - for (int32_t token : region_) { - // Ignore this token if it corresponds to an operator rather than a region. - if (token >= OP_UNION) - continue; - - // Calculate the distance to this surface. - // Note the off-by-one indexing - bool coincident {std::abs(token) == std::abs(on_surface)}; - double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)}; - - // Check if this distance is the new minimum. - if (d < min_dist) { - if (min_dist - d >= FP_PRECISION * min_dist) { - min_dist = d; - i_surf = -token; - } - } - } - - return {min_dist, i_surf}; -} - -//============================================================================== - void CSGCell::to_hdf5_inner(hid_t group_id) const { - write_string(group_id, "geom_type", "csg", false); - - // TODO: Move to Region::str() - // Write the region specification. - if (!region_.empty()) { - std::stringstream region_spec {}; - for (int32_t token : region_) { - if (token == OP_LEFT_PAREN) { - region_spec << " ("; - } else if (token == OP_RIGHT_PAREN) { - region_spec << " )"; - } else if (token == OP_COMPLEMENT) { - region_spec << " ~"; - } else if (token == OP_INTERSECTION) { - } else if (token == OP_UNION) { - region_spec << " |"; - } else { - // Note the off-by-one indexing - auto surf_id = model::surfaces[abs(token) - 1]->id_; - region_spec << " " << ((token > 0) ? surf_id : -surf_id); - } - } - write_string(group_id, "region", region_spec.str(), false); - } + write_string(group_id, "region", region_.str(), false); } -BoundingBox CSGCell::bounding_box_simple() const -{ - BoundingBox bbox; - for (int32_t token : region_) { - bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0); - } - return bbox; -} - -void CSGCell::apply_demorgan( - vector::iterator start, vector::iterator stop) -{ - do { - if (*start < OP_UNION) { - *start *= -1; - } else if (*start == OP_UNION) { - *start = OP_INTERSECTION; - } else if (*start == OP_INTERSECTION) { - *start = OP_UNION; - } - start++; - } while (start < stop); -} +//============================================================================== vector::iterator CSGCell::find_left_parenthesis( vector::iterator start, const vector& infix) @@ -779,75 +420,395 @@ vector::iterator CSGCell::find_left_parenthesis( return it; } -void CSGCell::remove_complement_ops(vector& infix) -{ - auto it = std::find(infix.begin(), infix.end(), OP_COMPLEMENT); - while (it != infix.end()) { - // Erase complement - infix.erase(it); +//============================================================================== +// Region implementation +//============================================================================== - // Define stop given left parenthesis or not - auto stop = it; - if (*it == OP_LEFT_PAREN) { - int depth = 1; - do { - stop++; - if (*stop > OP_COMPLEMENT) { - if (*stop == OP_RIGHT_PAREN) { - depth--; - } else { - depth++; - } +Region::Region(std::string region_spec, int32_t cell_id) +{ + // Check if region_spec is not empty. + if (!region_spec.empty()) { + // Parse all halfspaces and operators except for intersection (whitespace). + for (int i = 0; i < region_spec.size();) { + if (region_spec[i] == '(') { + expression_.push_back(OP_LEFT_PAREN); + i++; + + } else if (region_spec[i] == ')') { + expression_.push_back(OP_RIGHT_PAREN); + i++; + + } else if (region_spec[i] == '|') { + expression_.push_back(OP_UNION); + i++; + + } else if (region_spec[i] == '~') { + expression_.push_back(OP_COMPLEMENT); + i++; + + } else if (region_spec[i] == '-' || region_spec[i] == '+' || + std::isdigit(region_spec[i])) { + // This is the start of a halfspace specification. Iterate j until we + // find the end, then push-back everything between i and j. + int j = i + 1; + while (j < region_spec.size() && std::isdigit(region_spec[j])) { + j++; } - } while (depth > 0); - it++; + expression_.push_back(std::stoi(region_spec.substr(i, j - i))); + i = j; + + } else if (std::isspace(region_spec[i])) { + i++; + + } else { + auto err_msg = + fmt::format("Region specification contains invalid character, \"{}\"", + region_spec[i]); + fatal_error(err_msg); + } } - // apply DeMorgan's law to any surfaces/operators between these - // positions in the RPN - apply_demorgan(it, stop); - // update iterator position - it = std::find(infix.begin(), infix.end(), OP_COMPLEMENT); - } -} - -BoundingBox CSGCell::bounding_box_complex(vector postfix) -{ - vector stack(postfix.size()); - int i_stack = -1; - - for (auto& token : postfix) { - if (token == OP_UNION) { - stack[i_stack - 1] = stack[i_stack - 1] | stack[i_stack]; - i_stack--; - } else if (token == OP_INTERSECTION) { - stack[i_stack - 1] = stack[i_stack - 1] & stack[i_stack]; - i_stack--; - } else { - i_stack++; - stack[i_stack] = model::surfaces[abs(token) - 1]->bounding_box(token > 0); + // Add in intersection operators where a missing operator is needed. + int i = 0; + while (i < expression_.size() - 1) { + bool left_compat { + (expression_[i] < OP_UNION) || (expression_[i] == OP_RIGHT_PAREN)}; + bool right_compat {(expression_[i + 1] < OP_UNION) || + (expression_[i + 1] == OP_LEFT_PAREN) || + (expression_[i + 1] == OP_COMPLEMENT)}; + if (left_compat && right_compat) { + expression_.insert(expression_.begin() + i + 1, OP_INTERSECTION); + } + i++; } - } - Ensures(i_stack == 0); - return stack.front(); -} + // Remove complement operators using DeMorgan's laws + auto it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT); + while (it != expression_.end()) { + // Erase complement + expression_.erase(it); + + // Define stop given left parenthesis or not + auto stop = it; + if (*it == OP_LEFT_PAREN) { + int depth = 1; + do { + stop++; + if (*stop > OP_COMPLEMENT) { + if (*stop == OP_RIGHT_PAREN) { + depth--; + } else { + depth++; + } + } + } while (depth > 0); + it++; + } + + // apply DeMorgan's law to any surfaces/operators between these + // positions in the RPN + apply_demorgan(it, stop); + // update iterator position + it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT); + } + + // Convert user IDs to surface indices. + for (auto& r : expression_) { + if (r < OP_UNION) { + const auto& it {model::surface_map.find(abs(r))}; + if (it == model::surface_map.end()) { + throw std::runtime_error { + "Invalid surface ID " + std::to_string(abs(r)) + + " specified in region for cell " + std::to_string(cell_id) + "."}; + } + r = (r > 0) ? it->second + 1 : -(it->second + 1); + } + } + + // Check if this is a simple cell. + simple_ = true; + for (int32_t token : expression_) { + if (token == OP_UNION) { + simple_ = false; + // Ensure intersections have precedence over unions + add_precedence(); + break; + } + } + + // If this cell is simple, remove all the superfluous operator tokens. + if (simple_) { + for (auto it = expression_.begin(); it != expression_.end(); it++) { + if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) { + expression_.erase(it); + it--; + } + } + } + expression_.shrink_to_fit(); -BoundingBox CSGCell::bounding_box() const -{ - if (simple_) { - return bounding_box_simple(); } else { - auto postfix = generate_postfix(this->id_, this->region_); - return bounding_box_complex(postfix); + simple_ = true; } } //============================================================================== -bool CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const +void Region::apply_demorgan( + vector::iterator start, vector::iterator stop) { - for (int32_t token : region_) { + do { + if (*start < OP_UNION) { + *start *= -1; + } else if (*start == OP_UNION) { + *start = OP_INTERSECTION; + } else if (*start == OP_INTERSECTION) { + *start = OP_UNION; + } + start++; + } while (start < stop); +} + +//============================================================================== +//! Add precedence for infix regions so intersections have higher +//! precedence than unions using parentheses. +//============================================================================== + +std::vector::iterator Region::add_parentheses( + std::vector::iterator start) +{ + int32_t start_token = *start; + // Add left parenthesis + if (start_token == OP_INTERSECTION) { + start = expression_.insert(start - 1, OP_LEFT_PAREN); + } else { + start = expression_.insert(start + 1, OP_LEFT_PAREN); + } + start++; + + // Initialize return iterator + auto return_iterator = expression_.begin(); + + // Add right parenthesis + // While the start iterator is within the bounds of infix + while (start < expression_.end()) { + start++; + + // If the current token is an operator and is different than the start token + if (*start >= OP_UNION && *start != start_token) { + // Skip wrapped regions but save iterator position to check precedence and + // add right parenthesis, right parenthesis position depends on the + // operator, when the operator is a union then do not include the operator + // in the region, when the operator is an intersection then include the + // operato and next surface + if (*start == OP_LEFT_PAREN) { + return_iterator = start; + int depth = 1; + do { + start++; + if (*start > OP_COMPLEMENT) { + if (*start == OP_RIGHT_PAREN) { + depth--; + } else { + depth++; + } + } + } while (depth > 0); + } else { + start = expression_.insert( + start_token == OP_UNION ? start - 1 : start, OP_RIGHT_PAREN); + if (return_iterator == expression_.begin()) { + return_iterator = start - 1; + } + return return_iterator; + } + } + } + // If we get here a right parenthesis hasn't been placed, + // return iterator + expression_.push_back(OP_RIGHT_PAREN); + if (return_iterator == expression_.begin()) { + return_iterator = start - 1; + } + return return_iterator; +} + +//============================================================================== + +void Region::add_precedence() +{ + int32_t current_op = 0; + + for (auto it = expression_.begin(); it != expression_.end(); it++) { + int32_t token = *it; + + if (token == OP_UNION || token == OP_INTERSECTION) { + if (current_op == 0) { + // Set the current operator if is hasn't been set + current_op = token; + } else if (token != current_op) { + // If the current operator doesn't match the token, add parenthesis to + // assert precedence + it = add_parentheses(it); + current_op = 0; + } + } else if (token > OP_COMPLEMENT) { + // If the token is a parenthesis reset the current operator + current_op = 0; + } + } +} + +//============================================================================== +//! Convert infix region specification to Reverse Polish Notation (RPN) +//! +//! This function uses the shunting-yard algorithm. +//============================================================================== + +vector Region::generate_postfix(int32_t cell_id) const +{ + vector rpn; + vector stack; + + for (int32_t token : expression_) { + if (token < OP_UNION) { + // If token is not an operator, add it to output + rpn.push_back(token); + } else if (token < OP_RIGHT_PAREN) { + // Regular operators union, intersection, complement + while (stack.size() > 0) { + int32_t op = stack.back(); + + if (op < OP_RIGHT_PAREN && ((token == OP_COMPLEMENT && token < op) || + (token != OP_COMPLEMENT && token <= op))) { + // While there is an operator, op, on top of the stack, if the token + // is left-associative and its precedence is less than or equal to + // that of op or if the token is right-associative and its precedence + // is less than that of op, move op to the output queue and push the + // token on to the stack. Note that only complement is + // right-associative. + rpn.push_back(op); + stack.pop_back(); + } else { + break; + } + } + + stack.push_back(token); + + } else if (token == OP_LEFT_PAREN) { + // If the token is a left parenthesis, push it onto the stack + stack.push_back(token); + + } else { + // If the token is a right parenthesis, move operators from the stack to + // the output queue until reaching the left parenthesis. + for (auto it = stack.rbegin(); *it != OP_LEFT_PAREN; it++) { + // If we run out of operators without finding a left parenthesis, it + // means there are mismatched parentheses. + if (it == stack.rend()) { + fatal_error(fmt::format( + "Mismatched parentheses in region specification for cell {}", + cell_id)); + } + rpn.push_back(stack.back()); + stack.pop_back(); + } + + // Pop the left parenthesis. + stack.pop_back(); + } + } + + while (stack.size() > 0) { + int32_t op = stack.back(); + + // If the operator is a parenthesis it is mismatched. + if (op >= OP_RIGHT_PAREN) { + fatal_error(fmt::format( + "Mismatched parentheses in region specification for cell {}", cell_id)); + } + + rpn.push_back(stack.back()); + stack.pop_back(); + } + + return rpn; +} + +//============================================================================== + +std::string Region::str() const +{ + std::stringstream region_spec {}; + if (!expression_.empty()) { + for (int32_t token : expression_) { + if (token == OP_LEFT_PAREN) { + region_spec << " ("; + } else if (token == OP_RIGHT_PAREN) { + region_spec << " )"; + } else if (token == OP_COMPLEMENT) { + region_spec << " ~"; + } else if (token == OP_INTERSECTION) { + } else if (token == OP_UNION) { + region_spec << " |"; + } else { + // Note the off-by-one indexing + auto surf_id = model::surfaces[abs(token) - 1]->id_; + region_spec << " " << ((token > 0) ? surf_id : -surf_id); + } + } + } + return region_spec.str(); +} + +//============================================================================== + +std::pair Region::distance( + Position r, Direction u, int32_t on_surface, Particle* p) const +{ + double min_dist {INFTY}; + int32_t i_surf {std::numeric_limits::max()}; + + + for (int32_t token : expression_) { + // Ignore this token if it corresponds to an operator rather than a region. + if (token >= OP_UNION) + continue; + + // Calculate the distance to this surface. + // Note the off-by-one indexing + bool coincident {std::abs(token) == std::abs(on_surface)}; + double d {model::surfaces[abs(token) - 1]->distance(r, u, coincident)}; + + // Check if this distance is the new minimum. + if (d < min_dist) { + if (min_dist - d >= FP_PRECISION * min_dist) { + min_dist = d; + i_surf = -token; + } + } + } + + return {min_dist, i_surf}; +} + +//============================================================================== + +bool Region::contains(Position r, Direction u, int32_t on_surface) const +{ + if (simple_) { + return contains_simple(r, u, on_surface); + } else { + return contains_complex(r, u, on_surface); + } +} + +//============================================================================== + +bool Region::contains_simple(Position r, Direction u, int32_t on_surface) const +{ + for (int32_t token : expression_) { // Assume that no tokens are operators. Evaluate the sense of particle with // respect to the surface and see if the token matches the sense. If the // particle's surface attribute is set and matches the token, that @@ -868,14 +829,13 @@ bool CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const //============================================================================== -bool CSGCell::contains_complex( - Position r, Direction u, int32_t on_surface) const +bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const { bool in_cell = true; int total_depth = 0; // For each token - for (auto it = region_.begin(); it != region_.end(); it++) { + for (auto it = expression_.begin(); it != expression_.end(); it++) { int32_t token = *it; // If the token is a surface evaluate the sense @@ -926,6 +886,76 @@ bool CSGCell::contains_complex( return in_cell; } +//============================================================================== + +BoundingBox Region::bounding_box(int32_t cell_id) const +{ + if (simple_) { + return bounding_box_simple(); + } else { + auto postfix = generate_postfix(cell_id); + return bounding_box_complex(postfix); + } +} + +//============================================================================== + +BoundingBox Region::bounding_box_simple() const +{ + BoundingBox bbox; + for (int32_t token : expression_) { + bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0); + } + return bbox; +} + +//============================================================================== + +BoundingBox Region::bounding_box_complex(vector postfix) const +{ + vector stack(postfix.size()); + int i_stack = -1; + + for (auto& token : postfix) { + if (token == OP_UNION) { + stack[i_stack - 1] = stack[i_stack - 1] | stack[i_stack]; + i_stack--; + } else if (token == OP_INTERSECTION) { + stack[i_stack - 1] = stack[i_stack - 1] & stack[i_stack]; + i_stack--; + } else { + i_stack++; + stack[i_stack] = model::surfaces[abs(token) - 1]->bounding_box(token > 0); + } + } + + Ensures(i_stack == 0); + return stack.front(); +} + +//============================================================================== + +vector Region::surfaces() const +{ + if (simple_) { + return expression_; + } + + vector surfaces = expression_; + + auto it = std::find_if(surfaces.begin(), surfaces.end(), + [&](const auto& value) { return value >= OP_UNION; }); + + while (it != surfaces.end()) { + surfaces.erase(it); + + it = std::find_if(surfaces.begin(), surfaces.end(), + [&](const auto& value) { return value >= OP_UNION; }); + } + + return surfaces; +} + //============================================================================== // Non-method functions //============================================================================== diff --git a/src/geometry.cpp b/src/geometry.cpp index 216650ff18..29ca8b1bf8 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -425,7 +425,7 @@ BoundaryInfo distance_to_boundary(Particle& p) // positive half-space were given in the region specification. Thus, we // have to explicitly check which half-space the particle would be // traveling into if the surface is crossed - if (c.simple_) { + if (c.is_simple()) { info.surface_index = level_surf_cross; } else { Position r_hit = r + d_surf * u; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index c3458d4697..2614719840 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -154,9 +154,8 @@ void partition_universes() // Collect the set of surfaces in this universe. std::unordered_set surf_inds; for (auto i_cell : univ->cells_) { - for (auto token : model::cells[i_cell]->region_) { - if (token < OP_UNION) - surf_inds.insert(std::abs(token) - 1); + for (auto token : model::cells[i_cell]->surfaces()) { + surf_inds.insert(std::abs(token) - 1); } } diff --git a/src/universe.cpp b/src/universe.cpp index 873d03bb74..f8f9a82d89 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -98,13 +98,11 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) // 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]->region_) { - 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); - } + for (auto token : model::cells[i_cell]->surfaces()) { + 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); } } @@ -116,7 +114,7 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) 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_) { + if (!model::cells[i_cell]->is_simple()) { for (auto& p : partitions_) p.push_back(i_cell); continue; @@ -125,18 +123,16 @@ UniversePartitioner::UniversePartitioner(const Universe& univ) // 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]->region_) { - 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_; - } + for (auto token : model::cells[i_cell]->surfaces()) { + 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_; } } } From 169d6d93dcee45696bef96c4cb9b8ec8b5659791 Mon Sep 17 00:00:00 2001 From: myerspat Date: Wed, 5 Oct 2022 21:14:48 -0400 Subject: [PATCH 3/5] Added Doxygen comments to Region data and methods. Added default implementation for new void functions. --- include/openmc/cell.h | 58 +++++++++++++++++++++++++++++++++++++++---- 1 file changed, 53 insertions(+), 5 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 0891c062f3..248c2b23d7 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -62,34 +62,84 @@ public: //---------------------------------------------------------------------------- // Methods + + //! \brief Determine if a cell contains the particle at a given location. + //! + //! The bounds of the cell are determined by a logical expression involving + //! surface half-spaces. The expression used is given in infix notation + //! + //! The function is split into two cases, one for simple cells (those + //! involving only the intersection of half-spaces) and one for complex cells. + //! Both cases use short circuiting; however, in the case fo complex cells, + //! the complexity increases with the binary operators involved. + //! \param r The 3D Cartesian coordinate to check. + //! \param u A direction used to "break ties" the coordinates are very + //! close to a surface. + //! \param on_surface The signed index of a surface that the coordinate is + //! known to be on. This index takes precedence over surface sense + //! calculations. + bool contains(Position r, Direction u, int32_t on_surface) const; + + //! Find the oncoming boundary of this cell. std::pair distance( Position r, Direction u, int32_t on_surface, Particle* p) const; - bool contains(Position r, Direction u, int32_t on_surface) const; + + //! Get the BoundingBox for this cell. BoundingBox bounding_box(int32_t cell_id) const; + + //! Get the CSG expression as a string std::string str() const; + + //! Get a vector containing all the surfaces in the region expression vector surfaces() const; //---------------------------------------------------------------------------- // Accessors + + //! Get Boolean of if the cell is simple or not bool is_simple() const { return simple_; } private: //---------------------------------------------------------------------------- // Private Methods + + //! Get a vector of the region expression in postfix notation vector generate_postfix(int32_t cell_id) const; + + //! Determine if a particle is inside the cell for a simple cell (only + //! intersection operators) bool contains_simple(Position r, Direction u, int32_t on_surface) const; + + //! Determine if a particle is inside the cell for a complex cell. + //! + //! Uses the comobination of half-spaces and binary operators to determine + //! if short circuiting can be used. Short cicuiting uses the relative and + //! absolute depth of parenthases in the expression. bool contains_complex(Position r, Direction u, int32_t on_surface) const; + + //! BoundingBox if the paritcle is in a simple cell. BoundingBox bounding_box_simple() const; + + //! BoundingBox if the particle is in a complex cell. BoundingBox bounding_box_complex(vector postfix) const; + + //! Enfource precedence: Parenthases, Complement, Intersection, Union void add_precedence(); + + //! Add parenthesis to enforce precedence std::vector::iterator add_parentheses( std::vector::iterator start); + + //! Remove complement operators from the expression void remove_complement_ops(); + + //! Remove complement operators by using DeMorgan's laws void apply_demorgan( vector::iterator start, vector::iterator stop); //---------------------------------------------------------------------------- // Private Data + //! Definition of spatial region as Boolean expression of half-spaces // TODO: Should this be a vector of some other type vector expression_; @@ -98,8 +148,6 @@ private: //============================================================================== -// TODO: Think about what data members really need to live in this class versus -// putting them in the Region class class Cell { public: //---------------------------------------------------------------------------- @@ -155,10 +203,10 @@ public: virtual BoundingBox bounding_box() const = 0; //! Get a vector of surfaces in the cell - virtual vector surfaces() const = 0; + virtual vector surfaces() const { return vector(); } //! Check if the cell region expression is simple - virtual bool is_simple() const = 0; + virtual bool is_simple() const { return true; } //---------------------------------------------------------------------------- // Accessors From f60520e8c55e36c5dedb9d0ac52b1de035c02d0a Mon Sep 17 00:00:00 2001 From: myerspat Date: Thu, 6 Oct 2022 12:28:18 -0400 Subject: [PATCH 4/5] Removed simple_ from dagmc --- src/dagmc.cpp | 1 - 1 file changed, 1 deletion(-) diff --git a/src/dagmc.cpp b/src/dagmc.cpp index 0e04869fdc..a1b937e52c 100644 --- a/src/dagmc.cpp +++ b/src/dagmc.cpp @@ -557,7 +557,6 @@ DAGCell::DAGCell(std::shared_ptr dag_ptr, int32_t dag_idx) : Cell {}, dagmc_ptr_(dag_ptr), dag_index_(dag_idx) { geom_type_ = GeometryType::DAG; - simple_ = true; }; std::pair DAGCell::distance( From 647c9d185ea5d7a7be652da9fe22ebd5363f3ea5 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Fri, 4 Nov 2022 13:26:13 -0500 Subject: [PATCH 5/5] Update include/openmc/cell.h --- include/openmc/cell.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 248c2b23d7..39acaf0674 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -58,7 +58,7 @@ public: //---------------------------------------------------------------------------- // Constructors Region() {} - explicit Region(std::string region_expressioni, int32_t cell_id); + explicit Region(std::string region_spec, int32_t cell_id); //---------------------------------------------------------------------------- // Methods