From 630c14f12b3b4960cc0c1fe4a846c62102fc6ba6 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 2 Oct 2026 23:12:17 +0000 Subject: [PATCH 1/4] Represent cell regions as expression trees Regions were stored as infix token streams including parenthesis and operator tokens. Evaluating a complex region scanned the tokens while tracking parenthesis depth, precedence was enforced by inserting parentheses, and bounding boxes required converting the expression to postfix on every call. Region specifications are now parsed with a recursive descent parser into an expression tree of intersections and unions with half-spaces as leaves. Nested operators of the same type are merged, so the depth of the tree is the number of alternations between intersection and union. The tree is stored in pre-order with the index of the end of each subtree and of the parent of each node, so contains_complex evaluates it in a single loop without recursion, skipping the rest of a subtree as soon as its value is known. The half-spaces of a region are also stored as a contiguous list, which is all a simple region needs and is used for distance calculations. Complements are removed while parsing by applying De Morgan's laws to the tree. Previously complement operators were removed by flipping the operators in the token stream before precedence was enforced, so the complement of an expression mixing intersection and union without inner parentheses, such as ~(1 2 | 3), was interpreted as -1 | (-2 -3) instead of (-1 | -2) -3. Expressions written by the Python API are fully parenthesized and were not affected. A test comparing cell lookup against the Python region parser for random expressions is added. Region strings written to the summary file are generated from the tree. They are equivalent to the previous strings, but redundant parentheses are no longer kept. n_surfaces now returns the number of half-spaces rather than the number of tokens including operators. Results are identical for the tested models. Transport time is unchanged for lattice and simple geometries, and reduced by 16% for a water tank containing 400 rods, where the water region is the tank minus the rods. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014RN7JroEkrbkLKVHbMDap9 --- include/openmc/cell.h | 67 ++- src/cell.cpp | 672 ++++++++++-------------- tests/unit_tests/test_region_parsing.py | 78 +++ 3 files changed, 375 insertions(+), 442 deletions(-) create mode 100644 tests/unit_tests/test_region_parsing.py diff --git a/include/openmc/cell.h b/include/openmc/cell.h index c2a09d774b5..c0fd03d0053 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -66,12 +66,12 @@ class Region { //! \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 + //! surface half-spaces, stored as an expression tree of intersections and + //! unions with half-spaces as leaves. //! //! 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. + //! Both cases use short circuiting. //! \param r The 3D Cartesian coordinate to check. //! \param u A direction used to "break ties" the coordinates are very //! close to a surface. @@ -90,24 +90,36 @@ class Region { //! Get the CSG expression as a string std::string str() const; - //! Get a vector containing all the surfaces in the region expression + //! Get a vector containing all the half-spaces in the region expression vector surfaces() const; - //! Get size of surfaces - int n_surfaces() const { return expression_.size(); } + //! Get the number of half-spaces in the region expression + int n_surfaces() const; //---------------------------------------------------------------------------- // Accessors //! Get Boolean of if the cell is simple or not - bool is_simple() const { return simple_; } + bool is_simple() const { return nodes_.empty(); } private: //---------------------------------------------------------------------------- - // Private Methods + // Types + + //! Node of the region expression tree. Nodes are stored in pre-order, so + //! the children of an operator node follow it, and the subtree of a node + //! ends just before index end. Children of an operator node are never + //! operator nodes of the same type. + struct Node { + enum class Type : int8_t { HALFSPACE, INTERSECTION, UNION }; + Type type; + int32_t halfspace; //!< Signed surface index + 1 for HALFSPACE nodes + int32_t end; //!< Index one past the last node of the subtree + int32_t parent; //!< Index of the parent node (-1 for the root) + }; - //! Get a vector of the region expression in postfix notation - vector generate_postfix(int32_t cell_id) const; + //---------------------------------------------------------------------------- + // Private Methods //! Determine if a particle is inside the cell for a simple cell (only //! intersection operators) @@ -115,9 +127,8 @@ class Region { //! Determine if a particle is inside the cell for a complex cell. //! - //! Uses the combination of half-spaces and binary operators to determine - //! if short circuiting can be used. Short circuiting uses the relative and - //! absolute depth of parentheses in the expression. + //! Evaluates the expression tree, skipping the remaining children of an + //! operator node as soon as its value is known. bool contains_complex(Position r, Direction u, int32_t on_surface) const; //! Find the nearest intersection with any surface in the region expression. @@ -128,32 +139,16 @@ class Region { std::pair distance_complex( Position r, Direction u, int32_t on_surface) const; - //! BoundingBox if the particle 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; - - //! Enforce precedence between intersections and unions - void enforce_precedence(); - - //! Add parenthesis to enforce precedence - void add_parentheses(int64_t 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_; - bool simple_; //!< Does the region contain only intersections? + //! Signed surface indices + 1 of the half-spaces in the region expression, + //! in order. A simple region is the intersection of these half-spaces. + vector halfspaces_; + + //! Expression tree of a complex region in pre-order (empty for a simple + //! region) + vector nodes_; }; //============================================================================== diff --git a/src/cell.cpp b/src/cell.cpp index 8be8fb8d48f..4618c52d55a 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -579,26 +579,147 @@ vector::iterator CSGCell::find_left_parenthesis( // Region implementation //============================================================================== +namespace { + +//! Expression tree node used while parsing a region specification +struct ParseNode { + enum class Type { HALFSPACE, INTERSECTION, UNION }; + Type type; + int32_t halfspace {0}; + vector children; +}; + +//! Recursive descent parser for a tokenized region specification. +//! Intersection has higher precedence than union, and complement applies to the +//! half-space or parenthesized expression that follows it. Complements are +//! removed while parsing by applying De Morgan's laws, and nested operators of +//! the same type are merged. +class RegionParser { +public: + RegionParser(const vector& tokens, int32_t cell_id) + : tokens_(tokens), cell_id_(cell_id) + {} + + ParseNode parse() + { + ParseNode root = parse_union(false); + if (pos_ != tokens_.size()) { + if (tokens_[pos_] == OP_RIGHT_PAREN) + mismatched_parentheses(); + invalid(); + } + return root; + } + +private: + ParseNode parse_union(bool negate) + { + ParseNode node { + negate ? ParseNode::Type::INTERSECTION : ParseNode::Type::UNION}; + add_child(node, parse_intersection(negate)); + while (pos_ < tokens_.size() && tokens_[pos_] == OP_UNION) { + ++pos_; + add_child(node, parse_intersection(negate)); + } + return collapse(std::move(node)); + } + + ParseNode parse_intersection(bool negate) + { + ParseNode node { + negate ? ParseNode::Type::UNION : ParseNode::Type::INTERSECTION}; + add_child(node, parse_unary(negate)); + while (pos_ < tokens_.size() && tokens_[pos_] == OP_INTERSECTION) { + ++pos_; + add_child(node, parse_unary(negate)); + } + return collapse(std::move(node)); + } + + ParseNode parse_unary(bool negate) + { + if (pos_ >= tokens_.size()) + invalid(); + int32_t token = tokens_[pos_++]; + if (token == OP_COMPLEMENT) { + return parse_unary(!negate); + } else if (token == OP_LEFT_PAREN) { + ParseNode node = parse_union(negate); + if (pos_ >= tokens_.size() || tokens_[pos_] != OP_RIGHT_PAREN) + mismatched_parentheses(); + ++pos_; + return node; + } else if (token < OP_UNION) { + ParseNode node {ParseNode::Type::HALFSPACE}; + node.halfspace = negate ? -token : token; + return node; + } else if (token == OP_RIGHT_PAREN) { + mismatched_parentheses(); + } + invalid(); + } + + //! Add a child to an operator node, merging it into the node if it is an + //! operator node of the same type + static void add_child(ParseNode& parent, ParseNode child) + { + if (child.type == parent.type) { + for (auto& grandchild : child.children) + parent.children.push_back(std::move(grandchild)); + } else { + parent.children.push_back(std::move(child)); + } + } + + //! Replace an operator node with a single child by that child + static ParseNode collapse(ParseNode node) + { + if (node.children.size() == 1) + return std::move(node.children.front()); + return node; + } + + [[noreturn]] void mismatched_parentheses() const + { + fatal_error(fmt::format( + "Mismatched parentheses in region specification for cell {}", cell_id_)); + } + + [[noreturn]] void invalid() const + { + fatal_error( + fmt::format("Invalid region specification for cell {}", cell_id_)); + } + + const vector& tokens_; + int32_t cell_id_; + std::size_t pos_ {0}; +}; + +} // namespace + Region::Region(std::string region_spec, int32_t cell_id) { + vector tokens; + // 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); + tokens.push_back(OP_LEFT_PAREN); i++; } else if (region_spec[i] == ')') { - expression_.push_back(OP_RIGHT_PAREN); + tokens.push_back(OP_RIGHT_PAREN); i++; } else if (region_spec[i] == '|') { - expression_.push_back(OP_UNION); + tokens.push_back(OP_UNION); i++; } else if (region_spec[i] == '~') { - expression_.push_back(OP_COMPLEMENT); + tokens.push_back(OP_COMPLEMENT); i++; } else if (region_spec[i] == '-' || region_spec[i] == '+' || @@ -609,7 +730,7 @@ Region::Region(std::string region_spec, int32_t cell_id) while (j < region_spec.size() && std::isdigit(region_spec[j])) { j++; } - expression_.push_back(std::stoi(region_spec.substr(i, j - i))); + tokens.push_back(std::stoi(region_spec.substr(i, j - i))); i = j; } else if (std::isspace(region_spec[i])) { @@ -625,54 +746,20 @@ Region::Region(std::string region_spec, int32_t cell_id) // Add in intersection operators where a missing operator is needed. int i = 0; - while (i < expression_.size() - 1) { + while (i + 1 < tokens.size()) { 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)}; + (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) { - expression_.insert(expression_.begin() + i + 1, OP_INTERSECTION); + tokens.insert(tokens.begin() + i + 1, OP_INTERSECTION); } i++; } - // Remove complement operators using DeMorgan's laws - auto it = std::find(expression_.begin(), expression_.end(), OP_COMPLEMENT); - while (it != expression_.end()) { - // Erase complement. Note that erase invalidates the iterator, so we have - // to use the iterator it returns, which points to the token that - // followed the complement operator. - it = expression_.erase(it); - if (it == expression_.end()) - break; - - // 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_) { + for (auto& r : tokens) { if (r < OP_UNION) { const auto& it {model::surface_map.find(abs(r))}; if (it == model::surface_map.end()) { @@ -683,266 +770,81 @@ Region::Region(std::string region_spec, int32_t 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 - enforce_precedence(); - break; - } - } - - // If this cell is simple, remove all the superfluous operator tokens. - if (simple_) { - expression_.erase(std::remove_if(expression_.begin(), expression_.end(), - [](int32_t token) { - return token == OP_INTERSECTION || - token > OP_COMPLEMENT; - }), - expression_.end()); - } - expression_.shrink_to_fit(); - - } else { - simple_ = true; } -} - -//============================================================================== - -void Region::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); -} - -//============================================================================== -//! Add precedence for infix regions so intersections have higher -//! precedence than unions using parentheses. -//============================================================================== - -void Region::add_parentheses(int64_t start) -{ - int32_t start_token = expression_[start]; - // Add left parenthesis and set new position to be after parenthesis - if (start_token == OP_UNION) { - start += 2; - } - expression_.insert(expression_.begin() + start - 1, OP_LEFT_PAREN); - - // Add right parenthesis - // While the start iterator is within the bounds of infix - while (start + 1 < expression_.size()) { - start++; - - // If the current token is an operator and is different than the start token - if (expression_[start] >= OP_UNION && expression_[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 - // operator and next surface - if (expression_[start] == OP_LEFT_PAREN) { - int depth = 1; - do { - start++; - if (expression_[start] > OP_COMPLEMENT) { - if (expression_[start] == OP_RIGHT_PAREN) { - depth--; - } else { - depth++; - } - } - } while (depth > 0); - } else { - if (start_token == OP_UNION) { - --start; - } - expression_.insert(expression_.begin() + start, OP_RIGHT_PAREN); - return; - } - } - } - // If we get here a right parenthesis hasn't been placed - expression_.push_back(OP_RIGHT_PAREN); -} - -//============================================================================== -//! Add parentheses to enforce operator precedence in region expressions -//! -//! This function ensures that intersection operators have higher precedence -//! than union operators by adding parentheses where needed. For example: -//! "1 2 | 3" becomes "(1 2) | 3" -//! "1 | 2 3" becomes "1 | (2 3)" -//! -//! The algorithm uses stacks to track the current operator type and its -//! position at each parenthesis depth level. When it encounters a different -//! operator at the same depth, it adds parentheses to group the -//! higher-precedence operations. -//============================================================================== - -void Region::enforce_precedence() -{ - // Stack tracking the operator type at each depth (0 = no operator seen yet) - vector op_stack = {0}; - - // Stack tracking where the operator sequence started at each depth - vector pos_stack = {0}; - - for (int64_t i = 0; i < expression_.size(); ++i) { - int32_t token = expression_[i]; - - if (token == OP_LEFT_PAREN) { - // Entering a new parenthesis level - push new tracking state - op_stack.push_back(0); - pos_stack.push_back(0); - continue; - } else if (token == OP_RIGHT_PAREN) { - // Exiting a parenthesis level - pop tracking state (keep at least one) - if (op_stack.size() > 1) { - op_stack.pop_back(); - pos_stack.pop_back(); - } - continue; - } - if (token == OP_UNION || token == OP_INTERSECTION) { - if (op_stack.back() == 0) { - // First operator at this depth - record it and its position - op_stack.back() = token; - pos_stack.back() = i; - } else if (token != op_stack.back()) { - // Encountered a different operator at the same depth - need to add - // parentheses to enforce precedence. Intersection has higher - // precedence, so we parenthesize the intersection terms. - if (op_stack.back() == OP_INTERSECTION) { - add_parentheses(pos_stack.back()); - } else { - add_parentheses(i); - } + // An empty specification is the region containing all of space + if (tokens.empty()) + return; - // Restart the scan since we modified the expression - i = -1; // Will be incremented to 0 by the for loop - op_stack = {0}; - pos_stack = {0}; - } + // Parse the tokens into an expression tree and store it in pre-order + ParseNode root = RegionParser(tokens, cell_id).parse(); + bool simple = true; + auto append = [&](const ParseNode& node, int32_t parent, auto& self) -> void { + int32_t i = nodes_.size(); + Node::Type type; + switch (node.type) { + case ParseNode::Type::HALFSPACE: + type = Node::Type::HALFSPACE; + break; + case ParseNode::Type::INTERSECTION: + type = Node::Type::INTERSECTION; + break; + default: + type = Node::Type::UNION; + simple = false; } - } + nodes_.push_back({type, node.halfspace, 0, parent}); + for (const auto& child : node.children) + self(child, i, self); + nodes_[i].end = nodes_.size(); + }; + append(root, -1, append); + + // Store the half-spaces in order. A simple region is just their + // intersection, so its expression tree is not needed. + for (const auto& node : nodes_) { + if (node.type == Node::Type::HALFSPACE) + halfspaces_.push_back(node.halfspace); + } + if (simple) + nodes_.clear(); + nodes_.shrink_to_fit(); } -//============================================================================== -//! 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 +std::string Region::str() 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(); + std::string region_spec; + auto write = [&](int32_t i, bool parentheses, auto& self) -> void { + const Node& node = nodes_[i]; + if (node.type == Node::Type::HALFSPACE) { + // Note the off-by-one indexing + auto surf_id = model::surfaces[abs(node.halfspace) - 1]->id_; + region_spec += + fmt::format(" {}", (node.halfspace > 0) ? surf_id : -surf_id); + return; } - } - - 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)); + if (parentheses) + region_spec += " ("; + for (int32_t j = i + 1; j < node.end; j = nodes_[j].end) { + if (j > i + 1 && node.type == Node::Type::UNION) + region_spec += " |"; + self(j, true, self); } - - 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); - } + if (parentheses) + region_spec += " )"; + }; + if (!nodes_.empty()) { + write(0, false, write); + } else { + for (int32_t token : halfspaces_) { + // Note the off-by-one indexing + auto surf_id = model::surfaces[abs(token) - 1]->id_; + region_spec += fmt::format(" {}", (token > 0) ? surf_id : -surf_id); } } - return region_spec.str(); + return region_spec; } //============================================================================== @@ -950,7 +852,7 @@ std::string Region::str() const std::pair Region::distance( Position r, Direction u, int32_t on_surface) const { - if (simple_) { + if (nodes_.empty()) { return distance_to_nearest_surface(r, u, on_surface, false); } else { return distance_complex(r, u, on_surface); @@ -965,11 +867,7 @@ std::pair Region::distance_to_nearest_surface(Position r, 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; - + for (int32_t token : halfspaces_) { // Calculate the distance to this surface. // Note the off-by-one indexing bool coincident {std::abs(token) == std::abs(on_surface)}; @@ -1034,7 +932,7 @@ std::pair Region::distance_complex( bool Region::contains(Position r, Direction u, int32_t on_surface) const { - if (simple_) { + if (nodes_.empty()) { return contains_simple(r, u, on_surface); } else { return contains_complex(r, u, on_surface); @@ -1045,11 +943,11 @@ bool Region::contains(Position r, Direction u, int32_t on_surface) const 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 - // overrides the determination based on sense(). + for (int32_t token : halfspaces_) { + // 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 overrides the determination based on + // sense(). if (token == on_surface) { } else if (-token == on_surface) { return false; @@ -1068,129 +966,91 @@ bool Region::contains_simple(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 = expression_.begin(); it != expression_.end(); it++) { - int32_t token = *it; - - // If the token is a surface evaluate the sense - // If the token is a union or intersection check to - // short circuit - if (token < OP_UNION) { - if (token == on_surface) { - in_cell = true; - } else if (-token == on_surface) { - in_cell = false; + // Evaluate the expression tree without recursion: descend to the first + // half-space of each operator node and, after evaluating a half-space, move + // up the tree until reaching an operator node whose value is not yet known. + // The remaining children of an operator node are skipped as soon as one of + // them determines its value. + int32_t i = 0; + while (true) { + // Descend to the first half-space in this subtree + while (nodes_[i].type != Node::Type::HALFSPACE) + ++i; + + // Evaluate the half-space + bool value; + int32_t token = nodes_[i].halfspace; + if (token == on_surface) { + value = true; + } else if (-token == on_surface) { + value = false; + } else { + // Note the off-by-one indexing + value = model::surfaces[abs(token) - 1]->sense(r, u) == (token > 0); + } + + // Move up the tree until reaching a node with children left to evaluate + while (true) { + int32_t i_parent = nodes_[i].parent; + if (i_parent < 0) + return value; + const Node& parent = nodes_[i_parent]; + bool intersection = parent.type == Node::Type::INTERSECTION; + int32_t next = nodes_[i].end; + if (value != intersection || next == parent.end) { + // The value of the parent is known + i = i_parent; } else { - // Note the off-by-one indexing - bool sense = model::surfaces[abs(token) - 1]->sense(r, u); - in_cell = (sense == (token > 0)); - } - } else if ((token == OP_UNION && in_cell == true) || - (token == OP_INTERSECTION && in_cell == false)) { - // If the total depth is zero return - if (total_depth == 0) { - return in_cell; + i = next; + break; } - - total_depth--; - - // While the iterator is within the bounds of the vector - int depth = 1; - do { - // Get next token - it++; - int32_t next_token = *it; - - // If the token is an a parenthesis - if (next_token > OP_COMPLEMENT) { - // Adjust depth accordingly - if (next_token == OP_RIGHT_PAREN) { - depth--; - } else { - depth++; - } - } - } while (depth > 0); - } else if (token == OP_LEFT_PAREN) { - total_depth++; - } else if (token == OP_RIGHT_PAREN) { - total_depth--; } } - 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); + if (nodes_.empty()) { + BoundingBox bbox; + for (int32_t token : halfspaces_) { + bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0); + } + return bbox; } -} - -//============================================================================== -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; + auto box = [&](int32_t i, auto& self) -> BoundingBox { + const Node& node = nodes_[i]; + if (node.type == Node::Type::HALFSPACE) { + return model::surfaces[abs(node.halfspace) - 1]->bounding_box( + node.halfspace > 0); + } + BoundingBox bbox = self(i + 1, self); + for (int32_t j = nodes_[i + 1].end; j < node.end; j = nodes_[j].end) { + if (node.type == Node::Type::INTERSECTION) { + bbox &= self(j, self); + } else { + bbox = bbox | self(j, self); + } + } + return bbox; + }; + return box(0, box); } //============================================================================== -BoundingBox Region::bounding_box_complex(vector postfix) const +vector Region::surfaces() 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); - } - } - - assert(i_stack == 0); - return stack.front(); + return halfspaces_; } //============================================================================== -vector Region::surfaces() const +int Region::n_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; + return halfspaces_.size(); } //============================================================================== diff --git a/tests/unit_tests/test_region_parsing.py b/tests/unit_tests/test_region_parsing.py new file mode 100644 index 00000000000..89ba7aa915c --- /dev/null +++ b/tests/unit_tests/test_region_parsing.py @@ -0,0 +1,78 @@ +import random +import xml.etree.ElementTree as ET + +import numpy as np +import pytest +import openmc +import openmc.lib + + +SURFACE_IDS = [1, 2, 3, 4, 5, 6, 7] + + +def random_expression(rng, depth): + """Random region expression using minimal parentheses so that operator + precedence and complements of unparenthesized expressions are exercised""" + if depth == 0 or rng.random() < 0.25: + return f"{rng.choice(['', '-'])}{rng.choice(SURFACE_IDS)}" + op = rng.choice([' ', ' | ']) + terms = [random_expression(rng, depth - 1) for _ in range(rng.randint(2, 3))] + # A union only needs parentheses inside an intersection + if op == ' ': + terms = [f'({t})' if '|' in t else t for t in terms] + expr = op.join(terms) + if rng.random() < 0.3: + expr = f'~({expr})' + return expr + + +@pytest.mark.parametrize('seed', range(20)) +def test_region_expression_parsing(run_in_tmpdir, seed): + """Region expressions with minimal parentheses and complements are + interpreted the same way by the C++ and Python parsers""" + if seed == 0: + # Complement of an expression mixing intersection and union + expr = '~(1 3 | -5)' + else: + expr = random_expression(random.Random(seed), 4) + + surfaces = {s.id: s for s in [ + openmc.XPlane(-1.0, surface_id=1), + openmc.XPlane(2.0, surface_id=2), + openmc.YPlane(0.5, surface_id=3), + openmc.YPlane(-2.5, surface_id=4), + openmc.ZPlane(1.5, surface_id=5), + openmc.Sphere(x0=1.0, r=3.0, surface_id=6), + openmc.ZCylinder(y0=-1.0, r=2.0, surface_id=7), + openmc.Sphere(r=10.0, boundary_type='vacuum', surface_id=99), + ]} + region = openmc.Region.from_expression(expr, surfaces) + model = openmc.Model() + model.geometry = openmc.Geometry([ + openmc.Cell(cell_id=1, region=-surfaces[99] & region), + openmc.Cell(cell_id=2, region=-surfaces[99] & ~region), + ]) + model.settings.run_mode = 'fixed source' + model.settings.particles = 1 + model.settings.batches = 1 + model.export_to_model_xml() + + # Give OpenMC the raw expressions rather than the fully parenthesized + # strings written by the Python API + tree = ET.parse('model.xml') + for cell in tree.getroot().iter('cell'): + if cell.get('id') == '1': + cell.set('region', f'-99 {expr}') + else: + cell.set('region', f'-99 ~({expr})') + tree.write('model.xml') + + points = np.random.default_rng(seed).uniform(-5.7, 5.7, size=(300, 3)) + openmc.lib.init() + try: + for p in points: + cell, _ = openmc.lib.find_cell(p) + expected = 1 if tuple(p) in region else 2 + assert cell.id == expected, f'{expr} at {p}' + finally: + openmc.lib.finalize() From 80f73afecae88ee70cf7322af3879dbd53ec2617 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 2 Oct 2026 23:24:49 +0000 Subject: [PATCH 2/4] Remove code left unused by the region expression tree CSGCell::find_left_parenthesis searched token streams and had no callers, the cell ID argument of Region::bounding_box was only used for error messages when converting the expression to postfix, and the and includes are no longer used in cell.cpp. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014RN7JroEkrbkLKVHbMDap9 --- include/openmc/cell.h | 16 ++-------------- src/cell.cpp | 39 +-------------------------------------- 2 files changed, 3 insertions(+), 52 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index c0fd03d0053..94e923c4701 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -85,7 +85,7 @@ class Region { Position r, Direction u, int32_t on_surface) const; //! Get the BoundingBox for this cell. - BoundingBox bounding_box(int32_t cell_id) const; + BoundingBox bounding_box() const; //! Get the CSG expression as a string std::string str() const; @@ -442,10 +442,7 @@ class CSGCell : public Cell { return region_.contains(r, u, on_surface); } - BoundingBox bounding_box() const override - { - return region_.bounding_box(id_); - } + BoundingBox bounding_box() const override { return region_.bounding_box(); } void to_hdf5_inner(hid_t group_id) const override; @@ -453,15 +450,6 @@ class CSGCell : public Cell { virtual GeometryType geom_type() const override { return GeometryType::CSG; } -protected: - //! 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) - //! \param start Starting position of the search - //! \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 4618c52d55a..e07962b2b1b 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -6,8 +6,6 @@ #include #include #include -#include -#include #include #include @@ -540,41 +538,6 @@ void CSGCell::to_hdf5_inner(hid_t group_id) const write_string(group_id, "region", region_.str(), false); } -//============================================================================== - -vector::iterator CSGCell::find_left_parenthesis( - vector::iterator start, const vector& infix) -{ - // start search at zero - int parenthesis_level = 0; - auto it = start; - while (it != infix.begin()) { - // look at two tokens at a time - int32_t one = *it; - int32_t two = *(it - 1); - - // decrement parenthesis level if there are two adjacent surfaces - if (one < OP_UNION && two < OP_UNION) { - parenthesis_level--; - // increment if there are two adjacent operators - } else if (one >= OP_UNION && two >= OP_UNION) { - parenthesis_level++; - } - - // if the level gets to zero, return the position - if (parenthesis_level == 0) { - // move the iterator back one before leaving the loop - // so that all tokens in the parenthesis block are included - it--; - break; - } - - // continue loop, one token at a time - it--; - } - return it; -} - //============================================================================== // Region implementation //============================================================================== @@ -1010,7 +973,7 @@ bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const //============================================================================== -BoundingBox Region::bounding_box(int32_t cell_id) const +BoundingBox Region::bounding_box() const { if (nodes_.empty()) { BoundingBox bbox; From e38ad1bae446321ac08c80e19e38091a2cc41ff8 Mon Sep 17 00:00:00 2001 From: Claude Date: Fri, 2 Oct 2026 23:53:59 +0000 Subject: [PATCH 3/4] Update C++ region tests for expression tree output Region::str() is now generated from the expression tree, in which nested operators of the same type are merged and redundant parentheses are not kept, so the expected string for the issue #3685 case changes from " ( ( -1 2 ( -3 4 ) ) | ( -5 6 ) )" to the equivalent " ( -1 2 -3 4 ) | ( -5 6 )". Add cases checking that complements of mixed expressions are applied after grouping. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014RN7JroEkrbkLKVHbMDap9 --- tests/cpp_unit_tests/test_region.cpp | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/tests/cpp_unit_tests/test_region.cpp b/tests/cpp_unit_tests/test_region.cpp index ba784c842c7..b6c44d93995 100644 --- a/tests/cpp_unit_tests/test_region.cpp +++ b/tests/cpp_unit_tests/test_region.cpp @@ -108,9 +108,23 @@ TEST_CASE("Test region simplification") SECTION("Original bug case from issue #3685") { - // Input: "-1 2 (-3 4) | (-5 6)" was being incorrectly interpreted + // Input: "-1 2 (-3 4) | (-5 6)" was being incorrectly interpreted. + // Nested intersections are merged and redundant parentheses are dropped. auto region = openmc::Region("(-1 2 (-3 4) | (-5 6))", 0); - REQUIRE(region.str() == " ( ( -1 2 ( -3 4 ) ) | ( -5 6 ) )"); + REQUIRE(region.str() == " ( -1 2 -3 4 ) | ( -5 6 )"); + } + + SECTION("Complement of a mixed expression") + { + // The complement applies to the grouped expression (1 2) | 3 + auto region = openmc::Region("~(1 2 | 3)", 0); + REQUIRE(region.str() == " ( -1 | -2 ) -3"); + } + + SECTION("Complement of a parenthesized subexpression") + { + auto region = openmc::Region("4 ~(1 | 2 3)", 0); + REQUIRE(region.str() == " 4 -1 ( -2 | -3 )"); } SECTION("Simple union - no extra parentheses needed") From f3d451e9052f10e080c53be8d49f6d4a24b07f50 Mon Sep 17 00:00:00 2001 From: Claude Date: Sat, 3 Oct 2026 13:24:10 +0000 Subject: [PATCH 4/4] Keep complex-region data out of line in Region A simple region needs only its half-spaces, so the expression tree of a complex region is moved behind a pointer. This keeps Region, and the cells holding it, small for the common case of simple cells. Co-Authored-By: Claude Opus 5.5 Claude-Session: https://claude.ai/code/session_014RN7JroEkrbkLKVHbMDap9 --- include/openmc/cell.h | 13 ++++++++---- src/cell.cpp | 47 +++++++++++++++++++++++-------------------- 2 files changed, 34 insertions(+), 26 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 94e923c4701..39d7449e7db 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -100,7 +100,7 @@ class Region { // Accessors //! Get Boolean of if the cell is simple or not - bool is_simple() const { return nodes_.empty(); } + bool is_simple() const { return !complex_; } private: //---------------------------------------------------------------------------- @@ -146,9 +146,14 @@ class Region { //! in order. A simple region is the intersection of these half-spaces. vector halfspaces_; - //! Expression tree of a complex region in pre-order (empty for a simple - //! region) - vector nodes_; + //! Data needed only by complex regions, kept out of line so that regions, + //! and the cells holding them, stay small for simple cells + struct Complex { + vector nodes; //!< Expression tree in pre-order + }; + + //! Data of a complex region (null for a simple region) + unique_ptr complex_; }; //============================================================================== diff --git a/src/cell.cpp b/src/cell.cpp index e07962b2b1b..66c68784993 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -504,8 +504,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // Get a tokenized representation of the region specification and apply De // Morgans law - Region region(region_spec, id_); - region_ = region; + region_ = Region(region_spec, id_); // Read the translation vector. if (check_for_node(cell_node, "translation")) { @@ -663,6 +662,8 @@ class RegionParser { Region::Region(std::string region_spec, int32_t cell_id) { + vector nodes; + vector tokens; // Check if region_spec is not empty. @@ -743,7 +744,7 @@ Region::Region(std::string region_spec, int32_t cell_id) ParseNode root = RegionParser(tokens, cell_id).parse(); bool simple = true; auto append = [&](const ParseNode& node, int32_t parent, auto& self) -> void { - int32_t i = nodes_.size(); + int32_t i = nodes.size(); Node::Type type; switch (node.type) { case ParseNode::Type::HALFSPACE: @@ -756,22 +757,21 @@ Region::Region(std::string region_spec, int32_t cell_id) type = Node::Type::UNION; simple = false; } - nodes_.push_back({type, node.halfspace, 0, parent}); + nodes.push_back({type, node.halfspace, 0, parent}); for (const auto& child : node.children) self(child, i, self); - nodes_[i].end = nodes_.size(); + nodes[i].end = nodes.size(); }; append(root, -1, append); // Store the half-spaces in order. A simple region is just their // intersection, so its expression tree is not needed. - for (const auto& node : nodes_) { + for (const auto& node : nodes) { if (node.type == Node::Type::HALFSPACE) halfspaces_.push_back(node.halfspace); } - if (simple) - nodes_.clear(); - nodes_.shrink_to_fit(); + if (!simple) + complex_ = make_unique(Complex {std::move(nodes)}); } //============================================================================== @@ -780,7 +780,8 @@ std::string Region::str() const { std::string region_spec; auto write = [&](int32_t i, bool parentheses, auto& self) -> void { - const Node& node = nodes_[i]; + const auto& nodes = complex_->nodes; + const Node& node = nodes[i]; if (node.type == Node::Type::HALFSPACE) { // Note the off-by-one indexing auto surf_id = model::surfaces[abs(node.halfspace) - 1]->id_; @@ -790,7 +791,7 @@ std::string Region::str() const } if (parentheses) region_spec += " ("; - for (int32_t j = i + 1; j < node.end; j = nodes_[j].end) { + for (int32_t j = i + 1; j < node.end; j = nodes[j].end) { if (j > i + 1 && node.type == Node::Type::UNION) region_spec += " |"; self(j, true, self); @@ -798,7 +799,7 @@ std::string Region::str() const if (parentheses) region_spec += " )"; }; - if (!nodes_.empty()) { + if (complex_) { write(0, false, write); } else { for (int32_t token : halfspaces_) { @@ -815,7 +816,7 @@ std::string Region::str() const std::pair Region::distance( Position r, Direction u, int32_t on_surface) const { - if (nodes_.empty()) { + if (!complex_) { return distance_to_nearest_surface(r, u, on_surface, false); } else { return distance_complex(r, u, on_surface); @@ -895,7 +896,7 @@ std::pair Region::distance_complex( bool Region::contains(Position r, Direction u, int32_t on_surface) const { - if (nodes_.empty()) { + if (!complex_) { return contains_simple(r, u, on_surface); } else { return contains_complex(r, u, on_surface); @@ -934,15 +935,16 @@ bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const // up the tree until reaching an operator node whose value is not yet known. // The remaining children of an operator node are skipped as soon as one of // them determines its value. + const auto& nodes = complex_->nodes; int32_t i = 0; while (true) { // Descend to the first half-space in this subtree - while (nodes_[i].type != Node::Type::HALFSPACE) + while (nodes[i].type != Node::Type::HALFSPACE) ++i; // Evaluate the half-space bool value; - int32_t token = nodes_[i].halfspace; + int32_t token = nodes[i].halfspace; if (token == on_surface) { value = true; } else if (-token == on_surface) { @@ -954,12 +956,12 @@ bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const // Move up the tree until reaching a node with children left to evaluate while (true) { - int32_t i_parent = nodes_[i].parent; + int32_t i_parent = nodes[i].parent; if (i_parent < 0) return value; - const Node& parent = nodes_[i_parent]; + const Node& parent = nodes[i_parent]; bool intersection = parent.type == Node::Type::INTERSECTION; - int32_t next = nodes_[i].end; + int32_t next = nodes[i].end; if (value != intersection || next == parent.end) { // The value of the parent is known i = i_parent; @@ -975,7 +977,7 @@ bool Region::contains_complex(Position r, Direction u, int32_t on_surface) const BoundingBox Region::bounding_box() const { - if (nodes_.empty()) { + if (!complex_) { BoundingBox bbox; for (int32_t token : halfspaces_) { bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0); @@ -983,14 +985,15 @@ BoundingBox Region::bounding_box() const return bbox; } + const auto& nodes = complex_->nodes; auto box = [&](int32_t i, auto& self) -> BoundingBox { - const Node& node = nodes_[i]; + const Node& node = nodes[i]; if (node.type == Node::Type::HALFSPACE) { return model::surfaces[abs(node.halfspace) - 1]->bounding_box( node.halfspace > 0); } BoundingBox bbox = self(i + 1, self); - for (int32_t j = nodes_[i + 1].end; j < node.end; j = nodes_[j].end) { + for (int32_t j = nodes[i + 1].end; j < node.end; j = nodes[j].end) { if (node.type == Node::Type::INTERSECTION) { bbox &= self(j, self); } else {