From 29a741a441269d511b824e106b0b29a63493333a Mon Sep 17 00:00:00 2001 From: myerspat Date: Mon, 25 Jul 2022 11:58:16 -0500 Subject: [PATCH 01/17] added infix evaluation and removed conversion of infix to postfix --- src/cell.cpp | 79 +++++++++++++++++++++++++--------------------------- 1 file changed, 38 insertions(+), 41 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index f9786ce3eb..2e26e2c177 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -515,7 +515,8 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } // Convert the infix region spec to RPN. - rpn_ = generate_rpn(id_, region_); + //rpn_ = generate_rpn(id_, region_); + rpn_ = region_; // Check if this is a simple cell. simple_ = true; @@ -768,50 +769,46 @@ 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 { - // Make a stack of booleans. We don't know how big it needs to be, but we do - // know that rpn.size() is an upper-bound. - vector stack(rpn_.size()); - int i_stack = -1; + bool in_cell = true; + bool negate = false; + int paren_depth = 0; + int negate_depth = 0; for (int32_t token : rpn_) { - // If the token is a binary operator (intersection/union), apply it to - // the last two items on the stack. If the token is a unary operator - // (complement), apply it to the last item on the stack. - 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 if (token == OP_COMPLEMENT) { - stack[i_stack] = !stack[i_stack]; - } else { - // If the token is not an operator, 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(). - i_stack++; - if (token == on_surface) { - stack[i_stack] = true; - } else if (-token == on_surface) { - stack[i_stack] = false; - } else { - // Note the off-by-one indexing - bool sense = model::surfaces[abs(token) - 1]->sense(r, u); - stack[i_stack] = (sense == (token > 0)); - } - } - } + if (token < OP_UNION && in_cell == true) { + if (token == on_surface) { + in_cell = true; + } else if (-token == on_surface) { + in_cell = false; + } else { + // Note the off-by-one + bool sense = model::surfaces[abs(token) - 1]->sense(r, u); + in_cell = (sense == (token > 0)); + } + } else if (token == OP_UNION) { + if (in_cell == false) { + in_cell = true; + } else if (paren_depth == 0) { + break; + } + } else if (token == OP_RIGHT_PAREN) { + paren_depth--; + } else if (token == OP_LEFT_PAREN) { + paren_depth++; + } else if (token == OP_COMPLEMENT) { + negate = true; + negate_depth = paren_depth; + continue; + } else if (paren_depth == 0) { + break; + } - if (i_stack == 0) { - // The one remaining bool on the stack indicates whether the particle is - // in the cell. - return stack[i_stack]; - } else { - // This case occurs if there is no region specification since i_stack will - // still be -1. - return true; + if (negate == true && negate_depth == paren_depth) { + in_cell = !in_cell; + negate = false; + } } + return in_cell; } //============================================================================== From e4a10c8d05cb587a40a1711fed6cd55c4ddc756c Mon Sep 17 00:00:00 2001 From: myerspat Date: Tue, 26 Jul 2022 16:38:54 -0500 Subject: [PATCH 02/17] Added new algorithm for complex cell finding using infix --- src/cell.cpp | 120 +++++++++++++++++++++++++-------------------------- 1 file changed, 60 insertions(+), 60 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 2e26e2c177..578fb024da 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -514,14 +514,16 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } } - // Convert the infix region spec to RPN. - //rpn_ = generate_rpn(id_, region_); + // Convert the infix region spec to infix with no complements + // using De Morgan's law. + // TODO: Convert rpn to something that makes sense rpn_ = region_; + remove_complement_ops(rpn_); // Check if this is a simple cell. simple_ = true; for (int32_t token : rpn_) { - if ((token == OP_COMPLEMENT) || (token == OP_UNION)) { + if (token == OP_UNION) { simple_ = false; break; } @@ -529,16 +531,11 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // If this cell is simple, remove all the superfluous operator tokens. if (simple_) { - size_t i0 = 0; - size_t i1 = 0; - while (i1 < rpn_.size()) { - if (rpn_[i1] < OP_UNION) { - rpn_[i0] = rpn_[i1]; - ++i0; + for (auto it = rpn_.begin(); it != rpn_.end(); it++) { + if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) { + rpn_.erase(it); } - ++i1; } - rpn_.resize(i0); } rpn_.shrink_to_fit(); @@ -648,7 +645,7 @@ BoundingBox CSGCell::bounding_box_simple() const void CSGCell::apply_demorgan( vector::iterator start, vector::iterator stop) { - while (start < stop) { + do { if (*start < OP_UNION) { *start *= -1; } else if (*start == OP_UNION) { @@ -657,16 +654,16 @@ void CSGCell::apply_demorgan( *start = OP_UNION; } start++; - } + } while (start < stop); } vector::iterator CSGCell::find_left_parenthesis( - vector::iterator start, const vector& rpn) + vector::iterator start, const vector& infix) { // start search at zero int parenthesis_level = 0; auto it = start; - while (it != rpn.begin()) { + while (it != infix.begin()) { // look at two tokens at a time int32_t one = *it; int32_t two = *(it - 1); @@ -693,21 +690,25 @@ vector::iterator CSGCell::find_left_parenthesis( return it; } -void CSGCell::remove_complement_ops(vector& rpn) +void CSGCell::remove_complement_ops(vector& infix) { - 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); - vector tmp(left, it + 1); + auto it = std::find(infix.begin(), infix.end(), OP_COMPLEMENT); + while (it != infix.end()) { + // Erase complement + infix.erase(it); + + // Define stop given left parenthesis or not + auto stop = it; + if (*it == OP_LEFT_PAREN) { + stop = std::find(infix.begin(), infix.end(), OP_RIGHT_PAREN); + it++; + } // apply DeMorgan's law to any surfaces/operators between these // positions in the RPN - apply_demorgan(left, it); - // remove complement operator - rpn.erase(it); + apply_demorgan(it, stop); // update iterator position - it = std::find(rpn.begin(), rpn.end(), OP_COMPLEMENT); + it = std::find(infix.begin(), infix.end(), OP_COMPLEMENT); } } @@ -770,43 +771,42 @@ bool CSGCell::contains_complex( Position r, Direction u, int32_t on_surface) const { bool in_cell = true; - bool negate = false; - int paren_depth = 0; - int negate_depth = 0; + + // For each token + for (auto it = rpn_.begin(); it != rpn_.end(); it++) { + int32_t token = *it; - for (int32_t token : rpn_) { - if (token < OP_UNION && in_cell == true) { - if (token == on_surface) { - in_cell = true; - } else if (-token == on_surface) { - in_cell = false; - } else { - // Note the off-by-one - bool sense = model::surfaces[abs(token) - 1]->sense(r, u); - in_cell = (sense == (token > 0)); - } - } else if (token == OP_UNION) { - if (in_cell == false) { - in_cell = true; - } else if (paren_depth == 0) { - break; - } - } else if (token == OP_RIGHT_PAREN) { - paren_depth--; - } else if (token == OP_LEFT_PAREN) { - paren_depth++; - } else if (token == OP_COMPLEMENT) { - negate = true; - negate_depth = paren_depth; - continue; - } else if (paren_depth == 0) { - break; - } - - if (negate == true && negate_depth == paren_depth) { - in_cell = !in_cell; - negate = false; + // 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; + } 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)) { + // While the iterator is within the bounds of the vector + do { + // Get next token + it++; + int32_t next_token = *it; + + // If the next token is a left parenthesis skip until + // the next right parenthesis, if the token is a right + // parenthesis leave short circuiting + if (next_token == OP_LEFT_PAREN) { + it = std::find(it, rpn_.end(), OP_RIGHT_PAREN); + } else if (next_token == OP_RIGHT_PAREN) { + break; + } + } while (it < rpn_.end() - 1); + } } return in_cell; } From 67d01417b7c3092d7333702971da357b9d3c4e33 Mon Sep 17 00:00:00 2001 From: myerspat Date: Wed, 27 Jul 2022 11:18:02 -0500 Subject: [PATCH 03/17] Fixed errors due to iterators --- src/cell.cpp | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 578fb024da..edd4da84d9 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -700,7 +700,7 @@ void CSGCell::remove_complement_ops(vector& infix) // Define stop given left parenthesis or not auto stop = it; if (*it == OP_LEFT_PAREN) { - stop = std::find(infix.begin(), infix.end(), OP_RIGHT_PAREN); + stop = std::find(it, infix.end(), OP_RIGHT_PAREN); it++; } @@ -708,7 +708,7 @@ void CSGCell::remove_complement_ops(vector& infix) // positions in the RPN apply_demorgan(it, stop); // update iterator position - it = std::find(infix.begin(), infix.end(), OP_COMPLEMENT); + it = std::find(it, infix.end(), OP_COMPLEMENT); } } @@ -801,7 +801,7 @@ bool CSGCell::contains_complex( // the next right parenthesis, if the token is a right // parenthesis leave short circuiting if (next_token == OP_LEFT_PAREN) { - it = std::find(it, rpn_.end(), OP_RIGHT_PAREN); + it = std::find(it, rpn_.end() - 1, OP_RIGHT_PAREN); } else if (next_token == OP_RIGHT_PAREN) { break; } From 7dfb35a39a553bcb7abe06447a2d8e0a2ff6dc46 Mon Sep 17 00:00:00 2001 From: myerspat Date: Wed, 27 Jul 2022 12:15:46 -0500 Subject: [PATCH 04/17] removed rpn_ and changed to only region_ --- include/openmc/cell.h | 2 -- src/cell.cpp | 32 ++++++++++++++------------------ src/geometry_aux.cpp | 2 +- src/universe.cpp | 4 ++-- 4 files changed, 17 insertions(+), 23 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 01f0375865..8879017072 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -202,8 +202,6 @@ public: //! Definition of spatial region as Boolean expression of half-spaces vector region_; - //! Reverse Polish notation for region expression - vector rpn_; bool simple_; //!< Does the region contain only intersections? //! \brief Neighboring cells in the same universe. diff --git a/src/cell.cpp b/src/cell.cpp index edd4da84d9..943722f353 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -497,8 +497,10 @@ CSGCell::CSGCell(pugi::xml_node cell_node) region_spec = get_node_value(cell_node, "region"); } - // Get a tokenized representation of the region specification. + // Get a tokenized representation of the region specification + // and apply De Morgan's laws to remove complements. region_ = tokenize(region_spec); + remove_complement_ops(region_); region_.shrink_to_fit(); // Convert user IDs to surface indices. @@ -514,15 +516,9 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } } - // Convert the infix region spec to infix with no complements - // using De Morgan's law. - // TODO: Convert rpn to something that makes sense - rpn_ = region_; - remove_complement_ops(rpn_); - // Check if this is a simple cell. simple_ = true; - for (int32_t token : rpn_) { + for (int32_t token : region_) { if (token == OP_UNION) { simple_ = false; break; @@ -531,13 +527,13 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // If this cell is simple, remove all the superfluous operator tokens. if (simple_) { - for (auto it = rpn_.begin(); it != rpn_.end(); it++) { + for (auto it = region_.begin(); it != region_.end(); it++) { if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) { - rpn_.erase(it); + region_.erase(it); } } } - rpn_.shrink_to_fit(); + region_.shrink_to_fit(); // Read the translation vector. if (check_for_node(cell_node, "translation")) { @@ -581,7 +577,7 @@ std::pair CSGCell::distance( double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; - for (int32_t token : rpn_) { + for (int32_t token : region_) { // Ignore this token if it corresponds to an operator rather than a region. if (token >= OP_UNION) continue; @@ -636,7 +632,7 @@ void CSGCell::to_hdf5_inner(hid_t group_id) const BoundingBox CSGCell::bounding_box_simple() const { BoundingBox bbox; - for (int32_t token : rpn_) { + for (int32_t token : region_) { bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0); } return bbox; @@ -739,14 +735,14 @@ BoundingBox CSGCell::bounding_box_complex(vector rpn) BoundingBox CSGCell::bounding_box() const { - return simple_ ? bounding_box_simple() : bounding_box_complex(rpn_); + return simple_ ? bounding_box_simple() : bounding_box_complex(region_); } //============================================================================== bool CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const { - for (int32_t token : rpn_) { + for (int32_t token : region_) { // 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 @@ -773,7 +769,7 @@ bool CSGCell::contains_complex( bool in_cell = true; // For each token - for (auto it = rpn_.begin(); it != rpn_.end(); it++) { + for (auto it = region_.begin(); it != region_.end(); it++) { int32_t token = *it; // If the token is a surface evaluate the sense @@ -801,11 +797,11 @@ bool CSGCell::contains_complex( // the next right parenthesis, if the token is a right // parenthesis leave short circuiting if (next_token == OP_LEFT_PAREN) { - it = std::find(it, rpn_.end() - 1, OP_RIGHT_PAREN); + it = std::find(it, region_.end() - 1, OP_RIGHT_PAREN); } else if (next_token == OP_RIGHT_PAREN) { break; } - } while (it < rpn_.end() - 1); + } while (it < region_.end() - 1); } } return in_cell; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index ddce2469b2..c3458d4697 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -154,7 +154,7 @@ 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]->rpn_) { + for (auto token : model::cells[i_cell]->region_) { if (token < OP_UNION) surf_inds.insert(std::abs(token) - 1); } diff --git a/src/universe.cpp b/src/universe.cpp index f3f9cf8e8d..873d03bb74 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -98,7 +98,7 @@ 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]->rpn_) { + 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(); @@ -125,7 +125,7 @@ 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]->rpn_) { + 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)) { From 8a5df89ea5c899579deaa806f151519db2472aef Mon Sep 17 00:00:00 2001 From: myerspat Date: Thu, 28 Jul 2022 12:22:28 -0500 Subject: [PATCH 05/17] enforced precedence if intersections proceed unions without parenthesis --- include/openmc/cell.h | 2 + src/cell.cpp | 135 +++++++++++++++++++++++++++++++++--------- src/geometry_aux.cpp | 2 +- src/universe.cpp | 4 +- 4 files changed, 111 insertions(+), 32 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 8879017072..e812fc88f0 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -202,6 +202,8 @@ public: //! Definition of spatial region as Boolean expression of half-spaces vector region_; + //! Revised infix notation with no complements + vector region_no_complements_; bool simple_; //!< Does the region contain only intersections? //! \brief Neighboring cells in the same universe. diff --git a/src/cell.cpp b/src/cell.cpp index 943722f353..633cf1bc04 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -107,6 +107,66 @@ vector tokenize(const std::string region_spec) return tokens; } + +//============================================================================== +//! Add precedence for infix cell finding so intersections have higher +//! precedence than unions using parenthesis. +//============================================================================== + +std::vector::iterator +add_parenthesis(std::vector::iterator start, + std::vector &infix) { + // Add left parenthesis + start = infix.insert(start, OP_LEFT_PAREN); + start = start + 2; + + // Initialize return iterator + std::vector::iterator return_iterator = infix.end() - 1; + + // Add right parenthesis + // While the start iterator is within the bounds of infix + while (start < infix.end()) { + start++; + + // If we find a union or right parenthesis not wrapped by + // left and right parenthesis then place a right parenthesis, + // If we find a wrapped region return an iterator pointing to + // that wrapped region and continue looking to place a + // right parenthesis + if (*start == OP_UNION || *start == OP_RIGHT_PAREN) { + start = infix.insert(start, OP_RIGHT_PAREN); + return start - 1; + + } else if (*start == OP_LEFT_PAREN) { + return_iterator = start; + start = std::find(start, infix.end(), OP_RIGHT_PAREN); + } + } + + // If we get here a right parenthesis hasn't been placed, + // return iterator + infix.push_back(OP_RIGHT_PAREN); + 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 && current_op == 0) { + current_op = OP_UNION; + + } else if (current_op == OP_UNION && token == OP_INTERSECTION) { + it = add_parenthesis(it - 1, infix); + + } else if (token > OP_COMPLEMENT) { + current_op = 0; + } + } +} + //============================================================================== //! Convert infix region specification to Reverse Polish Notation (RPN) //! @@ -498,9 +558,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } // Get a tokenized representation of the region specification - // and apply De Morgan's laws to remove complements. region_ = tokenize(region_spec); - remove_complement_ops(region_); region_.shrink_to_fit(); // Convert user IDs to surface indices. @@ -516,9 +574,13 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } } + // Remove complement operators using De Morgan's law + region_no_complements_ = region_; + remove_complement_ops(region_no_complements_); + // Check if this is a simple cell. simple_ = true; - for (int32_t token : region_) { + for (int32_t token : region_no_complements_) { if (token == OP_UNION) { simple_ = false; break; @@ -527,13 +589,17 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // If this cell is simple, remove all the superfluous operator tokens. if (simple_) { - for (auto it = region_.begin(); it != region_.end(); it++) { + for (auto it = region_no_complements_.begin(); + it != region_no_complements_.end(); it++) { if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) { - region_.erase(it); + region_no_complements_.erase(it); } } + } else { + // Ensure intersections have precedence over unions + add_precedence(region_no_complements_); } - region_.shrink_to_fit(); + region_no_complements_.shrink_to_fit(); // Read the translation vector. if (check_for_node(cell_node, "translation")) { @@ -577,7 +643,7 @@ std::pair CSGCell::distance( double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; - for (int32_t token : region_) { + for (int32_t token : region_no_complements_) { // Ignore this token if it corresponds to an operator rather than a region. if (token >= OP_UNION) continue; @@ -632,7 +698,7 @@ void CSGCell::to_hdf5_inner(hid_t group_id) const BoundingBox CSGCell::bounding_box_simple() const { BoundingBox bbox; - for (int32_t token : region_) { + for (int32_t token : region_no_complements_) { bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0); } return bbox; @@ -735,14 +801,15 @@ BoundingBox CSGCell::bounding_box_complex(vector rpn) BoundingBox CSGCell::bounding_box() const { - return simple_ ? bounding_box_simple() : bounding_box_complex(region_); + return simple_ ? bounding_box_simple() + : bounding_box_complex(region_no_complements_); } //============================================================================== bool CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const { - for (int32_t token : region_) { + for (int32_t token : region_no_complements_) { // 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 @@ -767,13 +834,14 @@ bool CSGCell::contains_complex( Position r, Direction u, int32_t on_surface) const { bool in_cell = true; - + // For each token - for (auto it = region_.begin(); it != region_.end(); it++) { + for (auto it = region_no_complements_.begin(); + it != region_no_complements_.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 + // If the token is a union or intersection check to // short circuit if (token < OP_UNION) { if (token == on_surface) { @@ -794,14 +862,21 @@ bool CSGCell::contains_complex( int32_t next_token = *it; // If the next token is a left parenthesis skip until - // the next right parenthesis, if the token is a right - // parenthesis leave short circuiting - if (next_token == OP_LEFT_PAREN) { - it = std::find(it, region_.end() - 1, OP_RIGHT_PAREN); - } else if (next_token == OP_RIGHT_PAREN) { - break; + // the next right parenthesis, if the token is a right + // parenthesis leave short circuiting, if we go from + // intersections to union operators without parenthesis + // break shortc circuiting one behind the union + if (next_token >= OP_UNION) { + if (next_token == OP_LEFT_PAREN) { + it = std::find(it, region_no_complements_.end() - 1, OP_RIGHT_PAREN); + } else if (token == OP_RIGHT_PAREN) { + break; + } else if (token - next_token == 1) { + it--; + break; + } } - } while (it < region_.end() - 1); + } while (it < region_no_complements_.end() - 1); } } return in_cell; @@ -1108,7 +1183,8 @@ struct ParentCellStack { }; vector Cell::find_parent_cells( - int32_t instance, const Position& r) const { + int32_t instance, const Position& r) const +{ // create a temporary particle Particle dummy_particle {}; @@ -1118,8 +1194,8 @@ vector Cell::find_parent_cells( return find_parent_cells(instance, dummy_particle); } -vector Cell::find_parent_cells( - int32_t instance, Particle& p) const { +vector Cell::find_parent_cells(int32_t instance, Particle& p) const +{ // look up the particle's location exhaustive_find_cell(p); const auto& coords = p.coord(); @@ -1130,7 +1206,8 @@ vector Cell::find_parent_cells( for (auto it = coords.begin(); it != coords.end(); it++) { const auto& coord = *it; const auto& cell = model::cells[coord.cell]; - // if the cell at this level matches the current cell, stop adding to the stack + // if the cell at this level matches the current cell, stop adding to the + // stack if (coord.cell == model::cell_map[this->id_]) { cell_found = true; break; @@ -1141,7 +1218,8 @@ vector Cell::find_parent_cells( int lattice_idx = C_NONE; if (cell->type_ == Fill::LATTICE) { const auto& next_coord = *(it + 1); - lattice_idx = model::lattices[next_coord.lattice]->get_flat_index(next_coord.lattice_i); + lattice_idx = model::lattices[next_coord.lattice]->get_flat_index( + next_coord.lattice_i); } stack.push(coord.universe, {coord.cell, lattice_idx}); } @@ -1149,7 +1227,8 @@ vector Cell::find_parent_cells( // if this loop finished because the cell was found and // the instance matches the one requested in the call // we have the correct path and can return the stack - if (cell_found && stack.compute_instance(this->distribcell_index_) == instance) { + if (cell_found && + stack.compute_instance(this->distribcell_index_) == instance) { return stack.parent_cells(); } @@ -1157,9 +1236,7 @@ vector Cell::find_parent_cells( return exhaustive_find_parent_cells(instance); } - -vector Cell::exhaustive_find_parent_cells( - int32_t instance) const +vector Cell::exhaustive_find_parent_cells(int32_t instance) const { ParentCellStack stack; // start with this cell's universe diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index c3458d4697..881af15b76 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -154,7 +154,7 @@ 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_) { + for (auto token : model::cells[i_cell]->region_no_complements_) { if (token < OP_UNION) surf_inds.insert(std::abs(token) - 1); } diff --git a/src/universe.cpp b/src/universe.cpp index 873d03bb74..dc6f751c23 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -98,7 +98,7 @@ 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_) { + for (auto token : model::cells[i_cell]->region_no_complements_) { if (token < OP_UNION) { auto i_surf = std::abs(token) - 1; const auto* surf = model::surfaces[i_surf].get(); @@ -125,7 +125,7 @@ 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_) { + for (auto token : model::cells[i_cell]->region_no_complements_) { if (token < OP_UNION) { const auto* surf = model::surfaces[std::abs(token) - 1].get(); if (const auto* zplane = dynamic_cast(surf)) { From 9413e4699899dd9d19c6b583c67bcbfee2c86ea3 Mon Sep 17 00:00:00 2001 From: myerspat Date: Fri, 29 Jul 2022 11:15:53 -0500 Subject: [PATCH 06/17] changed break conditions for short ciruciting --- src/cell.cpp | 34 +++++++++++++++------------------- 1 file changed, 15 insertions(+), 19 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 633cf1bc04..f2e4444d4e 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -107,15 +107,14 @@ vector tokenize(const std::string region_spec) return tokens; } - //============================================================================== -//! Add precedence for infix cell finding so intersections have higher +//! Add precedence for infix cell finding so intersections have higher //! precedence than unions using parenthesis. //============================================================================== -std::vector::iterator -add_parenthesis(std::vector::iterator start, - std::vector &infix) { +std::vector::iterator add_parenthesis( + std::vector::iterator start, std::vector& infix) +{ // Add left parenthesis start = infix.insert(start, OP_LEFT_PAREN); start = start + 2; @@ -131,7 +130,7 @@ add_parenthesis(std::vector::iterator start, // If we find a union or right parenthesis not wrapped by // left and right parenthesis then place a right parenthesis, // If we find a wrapped region return an iterator pointing to - // that wrapped region and continue looking to place a + // that wrapped region and continue looking to place a // right parenthesis if (*start == OP_UNION || *start == OP_RIGHT_PAREN) { start = infix.insert(start, OP_RIGHT_PAREN); @@ -149,7 +148,8 @@ add_parenthesis(std::vector::iterator start, return return_iterator; } -void add_precedence(std::vector &infix) { +void add_precedence(std::vector& infix) +{ int32_t current_op = 0; for (auto it = infix.begin(); it != infix.end(); it++) { @@ -861,20 +861,16 @@ bool CSGCell::contains_complex( it++; int32_t next_token = *it; - // If the next token is a left parenthesis skip until - // the next right parenthesis, if the token is a right - // parenthesis leave short circuiting, if we go from - // intersections to union operators without parenthesis - // break shortc circuiting one behind the union - if (next_token >= OP_UNION) { + // If the next token is an operator and is not the same as the one that + // started the short circuiting, check if it is a left parenthesis to + // skip a section or break + if (next_token >= OP_UNION && token != next_token) { if (next_token == OP_LEFT_PAREN) { - it = std::find(it, region_no_complements_.end() - 1, OP_RIGHT_PAREN); - } else if (token == OP_RIGHT_PAREN) { + it = + std::find(it, region_no_complements_.end() - 1, OP_RIGHT_PAREN); + } else { break; - } else if (token - next_token == 1) { - it--; - break; - } + } } } while (it < region_no_complements_.end() - 1); } From d7904f7a137196b582b47b2e2d1b24f183ec7dd0 Mon Sep 17 00:00:00 2001 From: myerspat Date: Fri, 29 Jul 2022 13:29:21 -0500 Subject: [PATCH 07/17] now can handle multiple depths of parenthesis --- src/cell.cpp | 13 +++++++++++-- 1 file changed, 11 insertions(+), 2 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index f2e4444d4e..fc19c1c496 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -866,8 +866,17 @@ bool CSGCell::contains_complex( // skip a section or break if (next_token >= OP_UNION && token != next_token) { if (next_token == OP_LEFT_PAREN) { - it = - std::find(it, region_no_complements_.end() - 1, OP_RIGHT_PAREN); + int depth = 1; + while (depth > 0) { + it++; + if (*it > OP_COMPLEMENT) { + if (*it == OP_RIGHT_PAREN) { + depth--; + } else { + depth++; + } + } + } } else { break; } From bce3cc0da42651c01c7a63506f70b784769fff53 Mon Sep 17 00:00:00 2001 From: myerspat Date: Thu, 4 Aug 2022 15:46:10 -0400 Subject: [PATCH 08/17] Added rpn back for bounding box --- include/openmc/cell.h | 1 + src/cell.cpp | 22 +++++++++++++++------- 2 files changed, 16 insertions(+), 7 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index e812fc88f0..4d07b26199 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -204,6 +204,7 @@ public: vector region_; //! Revised infix notation with no complements vector region_no_complements_; + vector rpn_; bool simple_; //!< Does the region contain only intersections? //! \brief Neighboring cells in the same universe. diff --git a/src/cell.cpp b/src/cell.cpp index fc19c1c496..c89e437acd 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -138,7 +138,17 @@ std::vector::iterator add_parenthesis( } else if (*start == OP_LEFT_PAREN) { return_iterator = start; - start = std::find(start, infix.end(), OP_RIGHT_PAREN); + int depth = 1; + do { + start++; + if (*start > OP_COMPLEMENT) { + if (*start == OP_RIGHT_PAREN) { + depth--; + } else { + depth++; + } + } + } while (depth > 0); } } @@ -598,6 +608,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } else { // Ensure intersections have precedence over unions add_precedence(region_no_complements_); + rpn_ = generate_rpn(id_, region_no_complements_); } region_no_complements_.shrink_to_fit(); @@ -776,9 +787,6 @@ void CSGCell::remove_complement_ops(vector& infix) BoundingBox CSGCell::bounding_box_complex(vector rpn) { - // remove complements by adjusting surface signs and operators - remove_complement_ops(rpn); - vector stack(rpn.size()); int i_stack = -1; @@ -802,7 +810,7 @@ BoundingBox CSGCell::bounding_box_complex(vector rpn) BoundingBox CSGCell::bounding_box() const { return simple_ ? bounding_box_simple() - : bounding_box_complex(region_no_complements_); + : bounding_box_complex(rpn_); } //============================================================================== @@ -867,7 +875,7 @@ bool CSGCell::contains_complex( if (next_token >= OP_UNION && token != next_token) { if (next_token == OP_LEFT_PAREN) { int depth = 1; - while (depth > 0) { + do { it++; if (*it > OP_COMPLEMENT) { if (*it == OP_RIGHT_PAREN) { @@ -876,7 +884,7 @@ bool CSGCell::contains_complex( depth++; } } - } + } while (depth > 0); } else { break; } From a3c15fb07fcf51fec6fdeaa7faa665960b53c796 Mon Sep 17 00:00:00 2001 From: myerspat Date: Fri, 5 Aug 2022 17:01:09 -0400 Subject: [PATCH 09/17] Proper parenthesis handling for De Morgans --- src/cell.cpp | 12 +++++++++++- 1 file changed, 11 insertions(+), 1 deletion(-) diff --git a/src/cell.cpp b/src/cell.cpp index c89e437acd..12fb9f6740 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -773,7 +773,17 @@ void CSGCell::remove_complement_ops(vector& infix) // Define stop given left parenthesis or not auto stop = it; if (*it == OP_LEFT_PAREN) { - stop = std::find(it, infix.end(), OP_RIGHT_PAREN); + int depth = 1; + do { + stop++; + if (*stop > OP_COMPLEMENT) { + if (*stop == OP_RIGHT_PAREN) { + depth--; + } else { + depth++; + } + } + } while (depth > 0); it++; } From 809863c15ca760a708b62b5e844ebcb4aad2d3b7 Mon Sep 17 00:00:00 2001 From: myerspat Date: Mon, 15 Aug 2022 10:18:45 -0400 Subject: [PATCH 10/17] Removed second region expression so search finding relys on region_ --- include/openmc/cell.h | 9 ++------- src/cell.cpp | 45 +++++++++++++++++++++---------------------- src/geometry_aux.cpp | 2 +- src/universe.cpp | 4 ++-- 4 files changed, 27 insertions(+), 33 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 4d07b26199..0004a86c6f 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -160,14 +160,12 @@ protected: //! \param[in] instance of the cell to find parent cells for //! \param[in] p particle used to do a fast search for parent cells //! \return parent cells - vector find_parent_cells( - int32_t instance, Particle& p) const; + vector find_parent_cells(int32_t instance, Particle& p) const; //! Determine the path to this cell instance in the geometry hierarchy //! \param[in] instance of the cell to find parent cells for //! \return parent cells - vector exhaustive_find_parent_cells( - int32_t instance) const; + vector exhaustive_find_parent_cells(int32_t instance) const; //! Inner function for retrieving contained cells void get_contained_cells_inner( @@ -202,9 +200,6 @@ public: //! Definition of spatial region as Boolean expression of half-spaces vector region_; - //! Revised infix notation with no complements - vector region_no_complements_; - vector rpn_; bool simple_; //!< Does the region contain only intersections? //! \brief Neighboring cells in the same universe. diff --git a/src/cell.cpp b/src/cell.cpp index 12fb9f6740..8224ba2665 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -108,7 +108,7 @@ vector tokenize(const std::string region_spec) } //============================================================================== -//! Add precedence for infix cell finding so intersections have higher +//! Add precedence for infix regions so intersections have higher //! precedence than unions using parenthesis. //============================================================================== @@ -165,11 +165,17 @@ void add_precedence(std::vector& infix) for (auto it = infix.begin(); it != infix.end(); it++) { int32_t token = *it; + // If the token is a union another operator has not been found set the + // current operator to union + // If the current operator is a union and the token is an intersection + // assert precedence + // If the token is a parenthesis reset the current operator if (token == OP_UNION && current_op == 0) { current_op = OP_UNION; } else if (current_op == OP_UNION && token == OP_INTERSECTION) { it = add_parenthesis(it - 1, infix); + current_op = 0; } else if (token > OP_COMPLEMENT) { current_op = 0; @@ -567,8 +573,10 @@ CSGCell::CSGCell(pugi::xml_node cell_node) region_spec = get_node_value(cell_node, "region"); } - // Get a tokenized representation of the region specification + // Get a tokenized representation of the region specification and apply De + // Morgans law region_ = tokenize(region_spec); + remove_complement_ops(region_); region_.shrink_to_fit(); // Convert user IDs to surface indices. @@ -584,33 +592,26 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } } - // Remove complement operators using De Morgan's law - region_no_complements_ = region_; - remove_complement_ops(region_no_complements_); - // Check if this is a simple cell. simple_ = true; - for (int32_t token : region_no_complements_) { + for (int32_t token : region_) { if (token == OP_UNION) { simple_ = false; + // Ensure intersections have precedence over unions + add_precedence(region_); break; } } // If this cell is simple, remove all the superfluous operator tokens. if (simple_) { - for (auto it = region_no_complements_.begin(); - it != region_no_complements_.end(); it++) { + for (auto it = region_.begin(); it != region_.end(); it++) { if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) { - region_no_complements_.erase(it); + region_.erase(it); } } - } else { - // Ensure intersections have precedence over unions - add_precedence(region_no_complements_); - rpn_ = generate_rpn(id_, region_no_complements_); } - region_no_complements_.shrink_to_fit(); + region_.shrink_to_fit(); // Read the translation vector. if (check_for_node(cell_node, "translation")) { @@ -654,7 +655,7 @@ std::pair CSGCell::distance( double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; - for (int32_t token : region_no_complements_) { + for (int32_t token : region_) { // Ignore this token if it corresponds to an operator rather than a region. if (token >= OP_UNION) continue; @@ -709,7 +710,7 @@ void CSGCell::to_hdf5_inner(hid_t group_id) const BoundingBox CSGCell::bounding_box_simple() const { BoundingBox bbox; - for (int32_t token : region_no_complements_) { + for (int32_t token : region_) { bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0); } return bbox; @@ -819,15 +820,14 @@ BoundingBox CSGCell::bounding_box_complex(vector rpn) BoundingBox CSGCell::bounding_box() const { - return simple_ ? bounding_box_simple() - : bounding_box_complex(rpn_); + return simple_ ? bounding_box_simple() : bounding_box_complex(region_); } //============================================================================== bool CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const { - for (int32_t token : region_no_complements_) { + for (int32_t token : region_) { // 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 @@ -854,8 +854,7 @@ bool CSGCell::contains_complex( bool in_cell = true; // For each token - for (auto it = region_no_complements_.begin(); - it != region_no_complements_.end(); it++) { + for (auto it = region_.begin(); it != region_.end(); it++) { int32_t token = *it; // If the token is a surface evaluate the sense @@ -899,7 +898,7 @@ bool CSGCell::contains_complex( break; } } - } while (it < region_no_complements_.end() - 1); + } while (it < region_.end() - 1); } } return in_cell; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 881af15b76..c3458d4697 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -154,7 +154,7 @@ 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_no_complements_) { + for (auto token : model::cells[i_cell]->region_) { if (token < OP_UNION) surf_inds.insert(std::abs(token) - 1); } diff --git a/src/universe.cpp b/src/universe.cpp index dc6f751c23..873d03bb74 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -98,7 +98,7 @@ 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_no_complements_) { + 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(); @@ -125,7 +125,7 @@ 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_no_complements_) { + 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)) { From 82bbf7141352cac21295bac96d6e691692a170de Mon Sep 17 00:00:00 2001 From: myerspat Date: Mon, 15 Aug 2022 15:04:03 -0400 Subject: [PATCH 11/17] Reduced contains_complex to one while loop --- include/openmc/cell.h | 4 +- src/cell.cpp | 125 ++++++++++++++++++++++-------------------- src/geometry_aux.cpp | 2 +- src/universe.cpp | 4 +- 4 files changed, 72 insertions(+), 63 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 0004a86c6f..23cd0130a8 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -199,7 +199,9 @@ public: vector sqrtkT_; //! Definition of spatial region as Boolean expression of half-spaces - vector region_; + vector region_; + //! Region in infix that may be modified depending on simplicity + vector region_infix_; bool simple_; //!< Does the region contain only intersections? //! \brief Neighboring cells in the same universe. diff --git a/src/cell.cpp b/src/cell.cpp index 8224ba2665..a663d1144f 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -115,46 +115,62 @@ vector tokenize(const std::string region_spec) std::vector::iterator add_parenthesis( std::vector::iterator start, std::vector& infix) { + int32_t start_token = *start; // Add left parenthesis - start = infix.insert(start, OP_LEFT_PAREN); - start = start + 2; + 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 - std::vector::iterator return_iterator = infix.end() - 1; + std::vector::iterator return_iterator = infix.begin(); // Add right parenthesis // While the start iterator is within the bounds of infix while (start < infix.end()) { start++; - // If we find a union or right parenthesis not wrapped by - // left and right parenthesis then place a right parenthesis, - // If we find a wrapped region return an iterator pointing to - // that wrapped region and continue looking to place a - // right parenthesis - if (*start == OP_UNION || *start == OP_RIGHT_PAREN) { - start = infix.insert(start, OP_RIGHT_PAREN); - return start - 1; - - } else 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++; + // If the current token is an operator and is different than the start token + if (*start >= OP_UNION && *start != start_token) { + // Skip wraped regions but save iterator position to check precedence and + // add right parenthesis depending on the precedence of the original token + // when a right parenthesis is encountered of the opposite token + 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 if (start_token == OP_UNION) { + start = infix.insert(start - 1, OP_RIGHT_PAREN); + if (return_iterator == infix.begin()) { + return_iterator = start - 1; } - } while (depth > 0); + return return_iterator; + } else { + start = infix.insert(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; } @@ -168,15 +184,15 @@ void add_precedence(std::vector& infix) // If the token is a union another operator has not been found set the // current operator to union // If the current operator is a union and the token is an intersection - // assert precedence + // assert precedence // If the token is a parenthesis reset the current operator - if (token == OP_UNION && current_op == 0) { - current_op = OP_UNION; - - } else if (current_op == OP_UNION && token == OP_INTERSECTION) { - it = add_parenthesis(it - 1, infix); - current_op = 0; - + if (token == OP_UNION || token == OP_INTERSECTION) { + if (current_op == 0) { + current_op = token; + } else if (token != current_op) { + it = add_parenthesis(it, infix); + current_op = 0; + } } else if (token > OP_COMPLEMENT) { current_op = 0; } @@ -577,7 +593,6 @@ CSGCell::CSGCell(pugi::xml_node cell_node) // Morgans law region_ = tokenize(region_spec); remove_complement_ops(region_); - region_.shrink_to_fit(); // Convert user IDs to surface indices. for (auto& r : region_) { @@ -602,16 +617,18 @@ CSGCell::CSGCell(pugi::xml_node cell_node) break; } } + region_.shrink_to_fit(); + region_infix_ = region_; // If this cell is simple, remove all the superfluous operator tokens. if (simple_) { - for (auto it = region_.begin(); it != region_.end(); it++) { + for (auto it = region_infix_.begin(); it != region_infix_.end(); it++) { if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) { - region_.erase(it); + region_infix_.erase(it); } } + region_infix_.shrink_to_fit(); } - region_.shrink_to_fit(); // Read the translation vector. if (check_for_node(cell_node, "translation")) { @@ -655,7 +672,7 @@ std::pair CSGCell::distance( double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; - for (int32_t token : region_) { + for (int32_t token : region_infix_) { // Ignore this token if it corresponds to an operator rather than a region. if (token >= OP_UNION) continue; @@ -710,7 +727,7 @@ void CSGCell::to_hdf5_inner(hid_t group_id) const BoundingBox CSGCell::bounding_box_simple() const { BoundingBox bbox; - for (int32_t token : region_) { + for (int32_t token : region_infix_) { bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0); } return bbox; @@ -820,14 +837,14 @@ BoundingBox CSGCell::bounding_box_complex(vector rpn) BoundingBox CSGCell::bounding_box() const { - return simple_ ? bounding_box_simple() : bounding_box_complex(region_); + return simple_ ? bounding_box_simple() : bounding_box_complex(region_infix_); } //============================================================================== bool CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const { - for (int32_t token : region_) { + for (int32_t token : region_infix_) { // 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 @@ -854,7 +871,7 @@ bool CSGCell::contains_complex( bool in_cell = true; // For each token - for (auto it = region_.begin(); it != region_.end(); it++) { + for (auto it = region_infix_.begin(); it != region_infix_.end(); it++) { int32_t token = *it; // If the token is a surface evaluate the sense @@ -873,32 +890,22 @@ bool CSGCell::contains_complex( } else if ((token == OP_UNION && in_cell == true) || (token == OP_INTERSECTION && in_cell == false)) { // 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 next token is an operator and is not the same as the one that - // started the short circuiting, check if it is a left parenthesis to - // skip a section or break - if (next_token >= OP_UNION && token != next_token) { - if (next_token == OP_LEFT_PAREN) { - int depth = 1; - do { - it++; - if (*it > OP_COMPLEMENT) { - if (*it == OP_RIGHT_PAREN) { - depth--; - } else { - depth++; - } - } - } while (depth > 0); + // If the token is an a parenthesis + if (next_token > OP_COMPLEMENT) { + // Adjust depth accordingly + if (next_token == OP_RIGHT_PAREN) { + depth--; } else { - break; + depth++; } } - } while (it < region_.end() - 1); + } while (depth > 0 && it < region_infix_.end() - 1); } } return in_cell; diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index c3458d4697..7648d270ae 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -154,7 +154,7 @@ 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_) { + for (auto token : model::cells[i_cell]->region_infix_) { if (token < OP_UNION) surf_inds.insert(std::abs(token) - 1); } diff --git a/src/universe.cpp b/src/universe.cpp index 873d03bb74..723400ec05 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -98,7 +98,7 @@ 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_) { + for (auto token : model::cells[i_cell]->region_infix_) { if (token < OP_UNION) { auto i_surf = std::abs(token) - 1; const auto* surf = model::surfaces[i_surf].get(); @@ -125,7 +125,7 @@ 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_) { + for (auto token : model::cells[i_cell]->region_infix_) { if (token < OP_UNION) { const auto* surf = model::surfaces[std::abs(token) - 1].get(); if (const auto* zplane = dynamic_cast(surf)) { From 630117aab50172dd6659d36084b094dcf2148945 Mon Sep 17 00:00:00 2001 From: myerspat Date: Wed, 17 Aug 2022 08:19:22 -0400 Subject: [PATCH 12/17] Added return for a parenthesis depth of 0 --- src/cell.cpp | 14 +++++++++++++- 1 file changed, 13 insertions(+), 1 deletion(-) diff --git a/src/cell.cpp b/src/cell.cpp index a663d1144f..6b4fcafdaa 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -869,6 +869,7 @@ bool CSGCell::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_infix_.begin(); it != region_infix_.end(); it++) { @@ -889,6 +890,13 @@ bool CSGCell::contains_complex( } } 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; + } + + total_depth--; + // While the iterator is within the bounds of the vector int depth = 1; do { @@ -905,7 +913,11 @@ bool CSGCell::contains_complex( depth++; } } - } while (depth > 0 && it < region_infix_.end() - 1); + } while (depth > 0); + } else if (token == OP_LEFT_PAREN) { + total_depth++; + } else if (token == OP_RIGHT_PAREN) { + total_depth--; } } return in_cell; From 80c86925372bed9a59192b4e66b35d6917996550 Mon Sep 17 00:00:00 2001 From: myerspat Date: Thu, 25 Aug 2022 11:47:14 -0400 Subject: [PATCH 13/17] Removed one of the region vectors and changed HDF5 output for geometry --- include/openmc/cell.h | 2 -- src/cell.cpp | 49 +++++++++++++++++++++++-------------------- src/geometry_aux.cpp | 2 +- src/universe.cpp | 4 ++-- 4 files changed, 29 insertions(+), 28 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index 23cd0130a8..fd1883ce57 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -200,8 +200,6 @@ public: //! Definition of spatial region as Boolean expression of half-spaces vector region_; - //! Region in infix that may be modified depending on simplicity - vector region_infix_; bool simple_; //!< Does the region contain only intersections? //! \brief Neighboring cells in the same universe. diff --git a/src/cell.cpp b/src/cell.cpp index 6b4fcafdaa..4a61eb60e2 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -205,7 +205,7 @@ void add_precedence(std::vector& infix) //! This function uses the shunting-yard algorithm. //============================================================================== -vector generate_rpn(int32_t cell_id, vector infix) +vector generate_postfix(int32_t cell_id, vector infix) { vector rpn; vector stack; @@ -619,15 +619,14 @@ CSGCell::CSGCell(pugi::xml_node cell_node) } region_.shrink_to_fit(); - region_infix_ = region_; // If this cell is simple, remove all the superfluous operator tokens. if (simple_) { - for (auto it = region_infix_.begin(); it != region_infix_.end(); it++) { + for (auto it = region_.begin(); it != region_.end(); it++) { if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) { - region_infix_.erase(it); + region_.erase(it); } } - region_infix_.shrink_to_fit(); + region_.shrink_to_fit(); } // Read the translation vector. @@ -672,7 +671,7 @@ std::pair CSGCell::distance( double min_dist {INFTY}; int32_t i_surf {std::numeric_limits::max()}; - for (int32_t token : region_infix_) { + for (int32_t token : region_) { // Ignore this token if it corresponds to an operator rather than a region. if (token >= OP_UNION) continue; @@ -705,19 +704,18 @@ void CSGCell::to_hdf5_inner(hid_t group_id) const 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 { + if (token < OP_UNION) { // Note the off-by-one indexing auto surf_id = model::surfaces[abs(token) - 1]->id_; region_spec << " " << ((token > 0) ? surf_id : -surf_id); + } else if (token > OP_COMPLEMENT) { + if (token == OP_LEFT_PAREN) { + region_spec << " ("; + } else { + region_spec << " )"; + } + } else if (token == OP_UNION) { + region_spec << " |"; } } write_string(group_id, "region", region_spec.str(), false); @@ -727,7 +725,7 @@ void CSGCell::to_hdf5_inner(hid_t group_id) const BoundingBox CSGCell::bounding_box_simple() const { BoundingBox bbox; - for (int32_t token : region_infix_) { + for (int32_t token : region_) { bbox &= model::surfaces[abs(token) - 1]->bounding_box(token > 0); } return bbox; @@ -813,12 +811,12 @@ void CSGCell::remove_complement_ops(vector& infix) } } -BoundingBox CSGCell::bounding_box_complex(vector rpn) +BoundingBox CSGCell::bounding_box_complex(vector postfix) { - vector stack(rpn.size()); + vector stack(postfix.size()); int i_stack = -1; - for (auto& token : rpn) { + for (auto& token : postfix) { if (token == OP_UNION) { stack[i_stack - 1] = stack[i_stack - 1] | stack[i_stack]; i_stack--; @@ -837,14 +835,19 @@ BoundingBox CSGCell::bounding_box_complex(vector rpn) BoundingBox CSGCell::bounding_box() const { - return simple_ ? bounding_box_simple() : bounding_box_complex(region_infix_); + if (simple_) { + return bounding_box_simple(); + } else { + auto postfix = generate_postfix(this->id_, this->region_); + return bounding_box_complex(postfix); + } } //============================================================================== bool CSGCell::contains_simple(Position r, Direction u, int32_t on_surface) const { - for (int32_t token : region_infix_) { + for (int32_t token : region_) { // 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 @@ -872,7 +875,7 @@ bool CSGCell::contains_complex( int total_depth = 0; // For each token - for (auto it = region_infix_.begin(); it != region_infix_.end(); it++) { + for (auto it = region_.begin(); it != region_.end(); it++) { int32_t token = *it; // If the token is a surface evaluate the sense diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 7648d270ae..c3458d4697 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -154,7 +154,7 @@ 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_infix_) { + for (auto token : model::cells[i_cell]->region_) { if (token < OP_UNION) surf_inds.insert(std::abs(token) - 1); } diff --git a/src/universe.cpp b/src/universe.cpp index 723400ec05..873d03bb74 100644 --- a/src/universe.cpp +++ b/src/universe.cpp @@ -98,7 +98,7 @@ 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_infix_) { + 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(); @@ -125,7 +125,7 @@ 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_infix_) { + 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)) { From 090363b6372b32983fa2d4b23c84c3e4d2c83723 Mon Sep 17 00:00:00 2001 From: myerspat Date: Fri, 26 Aug 2022 12:47:08 -0400 Subject: [PATCH 14/17] Fixed error due to applying de morgans laws --- include/openmc/cell.h | 2 +- src/cell.cpp | 3 ++- 2 files changed, 3 insertions(+), 2 deletions(-) diff --git a/include/openmc/cell.h b/include/openmc/cell.h index fd1883ce57..8171967273 100644 --- a/include/openmc/cell.h +++ b/include/openmc/cell.h @@ -244,7 +244,7 @@ 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 rpn); + static BoundingBox bounding_box_complex(vector postfix); //! Applies DeMorgan's laws to a section of the RPN //! \param start Starting point for token modification diff --git a/src/cell.cpp b/src/cell.cpp index 4a61eb60e2..5c3e06eef3 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -624,6 +624,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) for (auto it = region_.begin(); it != region_.end(); it++) { if (*it == OP_INTERSECTION || *it > OP_COMPLEMENT) { region_.erase(it); + it--; } } region_.shrink_to_fit(); @@ -807,7 +808,7 @@ void CSGCell::remove_complement_ops(vector& infix) // positions in the RPN apply_demorgan(it, stop); // update iterator position - it = std::find(it, infix.end(), OP_COMPLEMENT); + it = std::find(infix.begin(), infix.end(), OP_COMPLEMENT); } } From bff8700c8b52dfbcf94f0aa65dc12cb53e56808d Mon Sep 17 00:00:00 2001 From: Patrick Myers <90068356+myerspat@users.noreply.github.com> Date: Tue, 30 Aug 2022 15:45:31 -0500 Subject: [PATCH 15/17] Simple wording and comment changes Co-authored-by: Paul Romano --- src/cell.cpp | 10 +++++----- 1 file changed, 5 insertions(+), 5 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 5c3e06eef3..9076b11364 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -109,10 +109,10 @@ vector tokenize(const std::string region_spec) //============================================================================== //! Add precedence for infix regions so intersections have higher -//! precedence than unions using parenthesis. +//! precedence than unions using parentheses. //============================================================================== -std::vector::iterator add_parenthesis( +std::vector::iterator add_parentheses( std::vector::iterator start, std::vector& infix) { int32_t start_token = *start; @@ -125,7 +125,7 @@ std::vector::iterator add_parenthesis( start++; // Initialize return iterator - std::vector::iterator return_iterator = infix.begin(); + auto return_iterator = infix.begin(); // Add right parenthesis // While the start iterator is within the bounds of infix @@ -134,7 +134,7 @@ std::vector::iterator add_parenthesis( // If the current token is an operator and is different than the start token if (*start >= OP_UNION && *start != start_token) { - // Skip wraped regions but save iterator position to check precedence and + // Skip wrapped regions but save iterator position to check precedence and // add right parenthesis depending on the precedence of the original token // when a right parenthesis is encountered of the opposite token if (*start == OP_LEFT_PAREN) { @@ -190,7 +190,7 @@ void add_precedence(std::vector& infix) if (current_op == 0) { current_op = token; } else if (token != current_op) { - it = add_parenthesis(it, infix); + it = add_parentheses(it, infix); current_op = 0; } } else if (token > OP_COMPLEMENT) { From 02254cc78787c3e21e3498c937a8fbdb2e7bc045 Mon Sep 17 00:00:00 2001 From: myerspat Date: Tue, 30 Aug 2022 17:00:59 -0400 Subject: [PATCH 16/17] Changed wording and rearranging --- src/cell.cpp | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index 9076b11364..bb524f0b3c 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -135,8 +135,10 @@ std::vector::iterator add_parentheses( // 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 depending on the precedence of the original token - // when a right parenthesis is encountered of the opposite token + // 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; @@ -150,14 +152,9 @@ std::vector::iterator add_parentheses( } } } while (depth > 0); - } else if (start_token == OP_UNION) { - start = infix.insert(start - 1, OP_RIGHT_PAREN); - if (return_iterator == infix.begin()) { - return_iterator = start - 1; - } - return return_iterator; } else { - start = infix.insert(start, OP_RIGHT_PAREN); + start = infix.insert( + start_token == OP_UNION ? start - 1 : start, OP_RIGHT_PAREN); if (return_iterator == infix.begin()) { return_iterator = start - 1; } @@ -181,19 +178,17 @@ void add_precedence(std::vector& infix) for (auto it = infix.begin(); it != infix.end(); it++) { int32_t token = *it; - // If the token is a union another operator has not been found set the - // current operator to union - // If the current operator is a union and the token is an intersection - // assert precedence - // If the token is a parenthesis reset the current operator 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; } } @@ -898,7 +893,7 @@ bool CSGCell::contains_complex( if (total_depth == 0) { return in_cell; } - + total_depth--; // While the iterator is within the bounds of the vector From 8650621f797f42d5482a4de4694273f907f344cf Mon Sep 17 00:00:00 2001 From: myerspat Date: Thu, 1 Sep 2022 17:01:01 -0400 Subject: [PATCH 17/17] Revereted hdf5 function back to original --- src/cell.cpp | 19 ++++++++++--------- 1 file changed, 10 insertions(+), 9 deletions(-) diff --git a/src/cell.cpp b/src/cell.cpp index bb524f0b3c..4ff35498c6 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -700,18 +700,19 @@ void CSGCell::to_hdf5_inner(hid_t group_id) const if (!region_.empty()) { std::stringstream region_spec {}; for (int32_t token : region_) { - if (token < OP_UNION) { + 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); - } else if (token > OP_COMPLEMENT) { - if (token == OP_LEFT_PAREN) { - region_spec << " ("; - } else { - region_spec << " )"; - } - } else if (token == OP_UNION) { - region_spec << " |"; } } write_string(group_id, "region", region_spec.str(), false);