diff --git a/include/openmc/cell.h b/include/openmc/cell.h index db8d60ed48..7206938c40 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -212,7 +212,26 @@ protected: bool contains_complex(Position r, Direction u, int32_t on_surface) const; BoundingBox bounding_box_simple() const; static BoundingBox bounding_box_complex(std::vector rpn); - static void apply_demorgan(std::vector& rpn); + + //! 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(std::vector::iterator start, + std::vector::iterator stop); + + //! Removes complement operators from the RPN + //! \param rpn The rpn to remove complement operators from. + static void remove_complement_ops(std::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) + //! \param start Starting position of the search + //! \param rpn The rpn being searched + static std::vector::iterator + find_left_parenthesis(std::vector::iterator start, + const std::vector& rpn); + }; //============================================================================== diff --git a/src/cell.cpp b/src/cell.cpp index cfb717c85e..3a7a2598fc 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -598,81 +598,89 @@ BoundingBox CSGCell::bounding_box_simple() const { return bbox; } +void CSGCell::apply_demorgan(std::vector::iterator start, + std::vector::iterator stop) +{ + while (start < stop) { + if (*start < OP_UNION) { *start *= -1; } + else if (*start == OP_UNION) { *start = OP_INTERSECTION; } + else if (*start == OP_INTERSECTION) { *start = OP_UNION; } + start++; + } +} -void CSGCell::apply_demorgan(std::vector& rpn) { - for (auto& token : rpn) { - if (token < OP_UNION) { token *= -1; } - else if (token == OP_UNION) { token = OP_INTERSECTION; } - else if (token == OP_INTERSECTION) { token = OP_UNION; } +std::vector::iterator +CSGCell::find_left_parenthesis(std::vector::iterator start, + const std::vector& rpn) { + // start search at zero + int parenthesis_level = 0; + auto it = start; + while (it != rpn.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; +} + +void CSGCell::remove_complement_ops(std::vector& rpn) { + auto it = std::find(rpn.begin(), rpn.end(), OP_COMPLEMENT); + while (it != rpn.end()) { + // find the opening parenthesis (if any) + auto left = find_left_parenthesis(it, rpn); + std::vector tmp(left, it+1); + + // apply DeMorgan's law to any surfaces/operators between these + // positions in the RPN + apply_demorgan(left, it); + // remove complement operator + rpn.erase(it); + // update iterator position + it = std::find(rpn.begin(), rpn.end(), OP_COMPLEMENT); } } BoundingBox CSGCell::bounding_box_complex(std::vector rpn) { + // remove complements by adjusting surface signs and operators + remove_complement_ops(rpn); - // if the last operator is a complement op, there is no - // sub-region that the complement connects to. This indicates - // that the entire region is a complement and we can apply - // De Morgan's laws immediately - if (rpn.back() == OP_COMPLEMENT) { - rpn.pop_back(); - apply_demorgan(rpn); - } + std::vector stack(rpn.size()); + int i_stack = -1; - // reverse the rpn to make popping easier - std::reverse(rpn.begin(), rpn.end()); - - BoundingBox current = model::surfaces[abs(rpn.back()) - 1]->bounding_box(rpn.back() > 0); - rpn.pop_back(); - - while (rpn.size()) { - // move through the rpn in twos - int32_t one = rpn.back(); rpn.pop_back(); - int32_t two = rpn.back(); rpn.pop_back(); - - // the first token should always be a surface - Expects(one < OP_UNION); - - if (two >= OP_UNION) { - if (two == OP_UNION) { - current |= model::surfaces[abs(one)-1]->bounding_box(one > 0); - } else if (two == OP_INTERSECTION) { - current &= model::surfaces[abs(one)-1]->bounding_box(one > 0); - } + for (auto& token : rpn) { + 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 { - // two surfaces in a row (left parenthesis), - // create sub-rpn for region in parenthesis - std::vector subrpn; - subrpn.push_back(one); - subrpn.push_back(two); - // add until last two tokens in the sub-rpn are operators - // (indicates a right parenthesis) - while (!((subrpn.back() >= OP_UNION) && (*(subrpn.rbegin() + 1) >= OP_UNION))) { - subrpn.push_back(rpn.back()); - rpn.pop_back(); - } - - // handle complement case using De Morgan's laws - if (subrpn.back() == OP_COMPLEMENT) { - subrpn.pop_back(); - apply_demorgan(subrpn); - subrpn.push_back(rpn.back()); - rpn.pop_back(); - } - // save the last operator, tells us how to combine this region - // with our current bounding box - int32_t op = subrpn.back(); subrpn.pop_back(); - // get bounding box for the subrpn - BoundingBox sub_box = bounding_box_complex(subrpn); - // combine the sub-rpn bounding box with our current cell box - if (op == OP_UNION) { - current |= sub_box; - } else if (op == OP_INTERSECTION) { - current &= sub_box; - } + i_stack++; + stack[i_stack] = model::surfaces[abs(token) - 1]->bounding_box(token > 0); } } - return current; + Ensures(i_stack == 0); + return stack.front(); } BoundingBox CSGCell::bounding_box() const { diff --git a/tests/unit_tests/test_complex_cell_bb.py b/tests/unit_tests/test_complex_cell_bb.py index 39431684d9..8db20f05c1 100644 --- a/tests/unit_tests/test_complex_cell_bb.py +++ b/tests/unit_tests/test_complex_cell_bb.py @@ -47,7 +47,7 @@ def complex_cell(run_in_tmpdir, mpi_intracomm): c1.region = ~(-s3 | +s4 | ~(+s13 & -s14)) c2 = openmc.Cell(fill=u238) - c2.region = +s2 & -s5 & +s12 & -s15 & ~(+s3 & -s4 & +s13 & -s14) + c2.region = ~(+s3 & -s4 & +s13 & -s14) & +s2 & -s5 & +s12 & -s15 c3 = openmc.Cell(fill=zr90) c3.region = ((+s1 & -s7 & +s17 & -s16) | (+s7 & -s6 & +s11 & -s17)) \