From 0f730d213ec983eef20e2e354fd3db1a85e527cb Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Sun, 8 Sep 2019 22:55:41 -0500
Subject: [PATCH 01/15] Added is_valid method to BoundingBox.

---
 include/openmc/surface.h | 7 +++++++
 1 file changed, 7 insertions(+)

diff --git a/include/openmc/surface.h b/include/openmc/surface.h
index d61ae07cef..9cdc460bb9 100644
--- a/include/openmc/surface.h
+++ b/include/openmc/surface.h
@@ -83,6 +83,13 @@ struct BoundingBox
     return *this;
   }
 
+  inline bool is_valid() {
+    return xmin <= xmax &&
+           ymin <= ymax &&
+           zmin <= zmax;
+  }
+
+
 };
 
 //==============================================================================

From 26fbbb6874d9f9cf11547c2551116d45e2f0c65c Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Sun, 8 Sep 2019 23:26:20 -0500
Subject: [PATCH 02/15] Moving to iterator-based method for advancing through
 the RPN.

---
 src/cell.cpp | 35 +++++++++++++++++++----------------
 1 file changed, 19 insertions(+), 16 deletions(-)

diff --git a/src/cell.cpp b/src/cell.cpp
index cfb717c85e..225d2a53df 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -618,16 +618,16 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
       apply_demorgan(rpn);
   }
 
-  // reverse the rpn to make popping easier
-  std::reverse(rpn.begin(), rpn.end());
+  auto it = rpn.begin();
+  BoundingBox current = model::surfaces[abs(*it) - 1]->bounding_box(*it > 0);
+  it++;
 
-  BoundingBox current = model::surfaces[abs(rpn.back()) - 1]->bounding_box(rpn.back() > 0);
-  rpn.pop_back();
-
-  while (rpn.size()) {
+  while (it < rpn.end()) {
     // move through the rpn in twos
-    int32_t one = rpn.back(); rpn.pop_back();
-    int32_t two = rpn.back(); rpn.pop_back();
+    int32_t one = *it;
+    it++;
+    int32_t two = *it;
+    it++;
 
     // the first token should always be a surface
     Expects(one < OP_UNION);
@@ -646,21 +646,24 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
       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();
+      while (!((*it >= OP_UNION) && (*(it + 1) >= OP_UNION))) {
+        subrpn.push_back(*it);
+        it++;
       }
 
+      // add first operator to the subrpn
+      subrpn.push_back(*it);
+      it++;
+
       // handle complement case using De Morgan's laws
-      if (subrpn.back() == OP_COMPLEMENT) {
-        subrpn.pop_back();
+      if (*it == OP_COMPLEMENT) {
         apply_demorgan(subrpn);
-        subrpn.push_back(rpn.back());
-        rpn.pop_back();
+        it++;
       }
       // save the last operator, tells us how to combine this region
       // with our current bounding box
-      int32_t op = subrpn.back(); subrpn.pop_back();
+      int32_t op = *it;
+      it++;
       // get bounding box for the subrpn
       BoundingBox sub_box = bounding_box_complex(subrpn);
       // combine the sub-rpn bounding box with our current cell box

From 85fa32150ad0788ed64ba1cd6011398d33bfd76f Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Mon, 9 Sep 2019 00:15:23 -0500
Subject: [PATCH 03/15] Still needs work

---
 include/openmc/cell.h |   3 +-
 src/cell.cpp          | 139 ++++++++++++++++++++++++++++++------------
 2 files changed, 101 insertions(+), 41 deletions(-)

diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index db8d60ed48..2c84955602 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -212,7 +212,8 @@ 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<int32_t> rpn);
-  static void apply_demorgan(std::vector<int32_t>& rpn);
+  static void apply_demorgan(std::vector<int32_t>::iterator start,
+                               std::vector<int32_t>::iterator end);
 };
 
 //==============================================================================
diff --git a/src/cell.cpp b/src/cell.cpp
index 225d2a53df..99ac9f43f4 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -598,12 +598,48 @@ BoundingBox CSGCell::bounding_box_simple() const {
   return bbox;
 }
 
+auto find_right_parenthesis(std::vector<int32_t>::iterator start,
+                            std::vector<int32_t>& rpn)
+{
+  // parenthesis level counter
+  int counter = 1;
+  std::cout << "FINDING" << std::endl;
+  // start at the beginning
+  auto it = start;
+  while(it < rpn.end()) {
+    int32_t one = *it;
+    int32_t two = *(it+1);
 
-void CSGCell::apply_demorgan(std::vector<int32_t>& rpn) {
-  for (auto& token : rpn) {
+    std::cout << "Counter: " << counter << std::endl;
+    std::cout << "One: " << one << std::endl;
+    std::cout << "Two: " << two << std::endl;
+
+    if (one < OP_UNION && two < OP_UNION) {
+      counter++;
+    } else if (one >= OP_UNION && two >= OP_UNION) {
+      counter--;
+    }
+
+    if (counter == 0) {
+      if (two == OP_COMPLEMENT) { it++; }
+      std::cout << "BREAK" << std::endl;
+      break;
+    }
+    it++;
+  }
+  std::cout << "DONE" << std::endl;
+  return it;
+}
+
+void CSGCell::apply_demorgan(std::vector<int32_t>::iterator start,
+                             std::vector<int32_t>::iterator end) {
+  auto it = start;
+  while (it <= end) {
+    int32_t token = *it;
     if (token < OP_UNION) { token *= -1; }
     else if (token == OP_UNION) { token = OP_INTERSECTION; }
     else if (token == OP_INTERSECTION) { token = OP_UNION; }
+    it++;
   }
 }
 
@@ -615,19 +651,26 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
   // De Morgan's laws immediately
   if (rpn.back() == OP_COMPLEMENT) {
       rpn.pop_back();
-      apply_demorgan(rpn);
+      apply_demorgan(rpn.begin(), rpn.end());
   }
 
-  auto it = rpn.begin();
-  BoundingBox current = model::surfaces[abs(*it) - 1]->bounding_box(*it > 0);
-  it++;
+  std::cout << "RPN" << std::endl;
+  for (auto i : rpn) { std::cout << i << " "; }
+  std::cout << std::endl;
 
+  // start with an invalid box
+  BoundingBox current;
+  current.xmin = INFTY;
+  current.xmax = -INFTY;
+
+  auto it = rpn.begin();
   while (it < rpn.end()) {
     // move through the rpn in twos
-    int32_t one = *it;
-    it++;
-    int32_t two = *it;
-    it++;
+    int32_t one = *it; it++;
+    int32_t two = *it; it++;
+
+    std::cout << "One: " << one << std::endl;
+    std::cout << "Two: " << two << std::endl;
 
     // the first token should always be a surface
     Expects(one < OP_UNION);
@@ -639,39 +682,55 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
         current &= model::surfaces[abs(one)-1]->bounding_box(one > 0);
       }
     } else {
-      // two surfaces in a row (left parenthesis),
-      // create sub-rpn for region in parenthesis
-      std::vector<int32_t> 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 (!((*it >= OP_UNION) && (*(it + 1) >= OP_UNION))) {
-        subrpn.push_back(*it);
-        it++;
-      }
 
-      // add first operator to the subrpn
-      subrpn.push_back(*it);
-      it++;
+      auto close = find_right_parenthesis(it, rpn);
 
-      // handle complement case using De Morgan's laws
-      if (*it == OP_COMPLEMENT) {
-        apply_demorgan(subrpn);
-        it++;
-      }
-      // save the last operator, tells us how to combine this region
-      // with our current bounding box
-      int32_t op = *it;
-      it++;
-      // 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;
+      std::vector<int32_t> subrpn(it - 2, close);
+
+      if (*close == OP_COMPLEMENT && *it < OP_UNION) {
+        apply_demorgan(it - 2, close);
+        it -= 2;
+      } else if (close == rpn.end()) {
+        current = model::surfaces[abs(one)-1]->bounding_box(one > 0);
+        it--;
+      } else {
+        int32_t op = *it;
+        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;
+        }
+        it = close;
       }
+      // // two surfaces in a row (left parenthesis),
+      // // create sub-rpn for region in parenthesis
+      // std::vector<int32_t> 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 (!((*it >= OP_UNION) && (*(it + 1) >= OP_UNION))) {
+      //   subrpn.push_back(*it);
+      //   it++;
+      // }
+
+      // // add first operator to the subrpn
+      // subrpn.push_back(*it);
+      // it++;
+
+      // // handle complement case using De Morgan's laws
+      // if (*it == OP_COMPLEMENT) {
+      //   apply_demorgan(subrpn);
+      //   it++;
+      // }
+      // // save the last operator, tells us how to combine this region
+      // // with our current bounding box
+      // int32_t op = *it;
+      // it++;
+      // // get bounding box for the subrpn
+      // BoundingBox sub_box = bounding_box_complex(subrpn);
     }
   }
 

From 0d02f8b2d29c0a343eff560f2492ab23f7d13ba0 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Mon, 9 Sep 2019 16:11:48 -0500
Subject: [PATCH 04/15] Working now. Needs cleanup.

---
 include/openmc/cell.h                    |   2 +-
 src/cell.cpp                             | 144 ++++++++++-------------
 tests/unit_tests/test_complex_cell_bb.py |   2 +-
 3 files changed, 62 insertions(+), 86 deletions(-)

diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 2c84955602..471b5c4f85 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -213,7 +213,7 @@ protected:
   BoundingBox bounding_box_simple() const;
   static BoundingBox bounding_box_complex(std::vector<int32_t> rpn);
   static void apply_demorgan(std::vector<int32_t>::iterator start,
-                               std::vector<int32_t>::iterator end);
+                             std::vector<int32_t>::iterator stop);
 };
 
 //==============================================================================
diff --git a/src/cell.cpp b/src/cell.cpp
index 99ac9f43f4..40615c1e89 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -598,53 +598,45 @@ BoundingBox CSGCell::bounding_box_simple() const {
   return bbox;
 }
 
+
+void CSGCell::apply_demorgan(std::vector<int32_t>::iterator start,
+                             std::vector<int32_t>::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++;
+  }
+}
+
 auto find_right_parenthesis(std::vector<int32_t>::iterator start,
                             std::vector<int32_t>& rpn)
 {
   // parenthesis level counter
-  int counter = 1;
-  std::cout << "FINDING" << std::endl;
+  int counter = 0;
+
   // start at the beginning
   auto it = start;
   while(it < rpn.end()) {
     int32_t one = *it;
     int32_t two = *(it+1);
 
-    std::cout << "Counter: " << counter << std::endl;
-    std::cout << "One: " << one << std::endl;
-    std::cout << "Two: " << two << std::endl;
-
     if (one < OP_UNION && two < OP_UNION) {
       counter++;
     } else if (one >= OP_UNION && two >= OP_UNION) {
       counter--;
+      while(*it >= OP_UNION) { it++; }
+      it--;
     }
 
-    if (counter == 0) {
-      if (two == OP_COMPLEMENT) { it++; }
-      std::cout << "BREAK" << std::endl;
-      break;
-    }
+    if (counter == 0) { break; }
     it++;
   }
-  std::cout << "DONE" << std::endl;
   return it;
 }
 
-void CSGCell::apply_demorgan(std::vector<int32_t>::iterator start,
-                             std::vector<int32_t>::iterator end) {
-  auto it = start;
-  while (it <= end) {
-    int32_t token = *it;
-    if (token < OP_UNION) { token *= -1; }
-    else if (token == OP_UNION) { token = OP_INTERSECTION; }
-    else if (token == OP_INTERSECTION) { token = OP_UNION; }
-    it++;
-  }
-}
-
 BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> 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
@@ -654,23 +646,23 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
       apply_demorgan(rpn.begin(), rpn.end());
   }
 
-  std::cout << "RPN" << std::endl;
-  for (auto i : rpn) { std::cout << i << " "; }
-  std::cout << std::endl;
-
-  // start with an invalid box
-  BoundingBox current;
-  current.xmin = INFTY;
-  current.xmax = -INFTY;
+  // check for parenthesis with complement at the start of the rpn
+  auto initial_right = find_right_parenthesis(rpn.begin(), rpn);
+  if (initial_right != rpn.end() && *initial_right == OP_COMPLEMENT) {
+    apply_demorgan(rpn.begin(), initial_right);
+    rpn.erase(initial_right);
+  }
 
   auto it = rpn.begin();
+  BoundingBox current = model::surfaces[abs(*it) - 1]->bounding_box(*it > 0);
+  it++;
+
   while (it < rpn.end()) {
     // move through the rpn in twos
-    int32_t one = *it; it++;
-    int32_t two = *it; it++;
-
-    std::cout << "One: " << one << std::endl;
-    std::cout << "Two: " << two << std::endl;
+    int32_t one = *it;
+    it++;
+    int32_t two = *it;
+    it++;
 
     // the first token should always be a surface
     Expects(one < OP_UNION);
@@ -682,55 +674,39 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
         current &= model::surfaces[abs(one)-1]->bounding_box(one > 0);
       }
     } else {
-
-      auto close = find_right_parenthesis(it, rpn);
-
-      std::vector<int32_t> subrpn(it - 2, close);
-
-      if (*close == OP_COMPLEMENT && *it < OP_UNION) {
-        apply_demorgan(it - 2, close);
-        it -= 2;
-      } else if (close == rpn.end()) {
-        current = model::surfaces[abs(one)-1]->bounding_box(one > 0);
-        it--;
-      } else {
-        int32_t op = *it;
-        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;
-        }
-        it = close;
+      // two surfaces in a row (left parenthesis),
+      // create sub-rpn for region in parenthesis
+      std::vector<int32_t> 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 (!((*it >= OP_UNION) && (*(it + 1) >= OP_UNION))) {
+        subrpn.push_back(*it);
+        it++;
       }
-      // // two surfaces in a row (left parenthesis),
-      // // create sub-rpn for region in parenthesis
-      // std::vector<int32_t> 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 (!((*it >= OP_UNION) && (*(it + 1) >= OP_UNION))) {
-      //   subrpn.push_back(*it);
-      //   it++;
-      // }
 
-      // // add first operator to the subrpn
-      // subrpn.push_back(*it);
-      // it++;
+      // add first operator to the subrpn
+      subrpn.push_back(*it);
+      it++;
 
-      // // handle complement case using De Morgan's laws
-      // if (*it == OP_COMPLEMENT) {
-      //   apply_demorgan(subrpn);
-      //   it++;
-      // }
-      // // save the last operator, tells us how to combine this region
-      // // with our current bounding box
-      // int32_t op = *it;
-      // it++;
-      // // get bounding box for the subrpn
-      // BoundingBox sub_box = bounding_box_complex(subrpn);
+      // handle complement case using De Morgan's laws
+      if (*it == OP_COMPLEMENT) {
+        apply_demorgan(subrpn.begin(), subrpn.end());
+        it++;
+      }
+      // save the last operator, tells us how to combine this region
+      // with our current bounding box
+      int32_t op = *it;
+      it++;
+      // 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;
+      }
     }
   }
 
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)) \

From 789d55a35a6ef75e3aa67aac76bd3cbed60e24d7 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Wed, 11 Sep 2019 09:32:26 -0500
Subject: [PATCH 05/15] Removing complement operator in a separate step.

---
 include/openmc/cell.h |  4 +++
 src/cell.cpp          | 84 ++++++++++++++++++++++---------------------
 2 files changed, 48 insertions(+), 40 deletions(-)

diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 471b5c4f85..90c03e3a90 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -214,6 +214,10 @@ protected:
   static BoundingBox bounding_box_complex(std::vector<int32_t> rpn);
   static void apply_demorgan(std::vector<int32_t>::iterator start,
                              std::vector<int32_t>::iterator stop);
+  static void remove_complements(std::vector<int32_t>& rpn);
+  static std::vector<int32_t>::iterator
+  find_left_parenthesis(std::vector<int32_t>::iterator start,
+                        const std::vector<int32_t>& rpn);
 };
 
 //==============================================================================
diff --git a/src/cell.cpp b/src/cell.cpp
index 40615c1e89..f33926b303 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -598,7 +598,6 @@ BoundingBox CSGCell::bounding_box_simple() const {
   return bbox;
 }
 
-
 void CSGCell::apply_demorgan(std::vector<int32_t>::iterator start,
                              std::vector<int32_t>::iterator stop)
 {
@@ -610,48 +609,62 @@ void CSGCell::apply_demorgan(std::vector<int32_t>::iterator start,
   }
 }
 
-auto find_right_parenthesis(std::vector<int32_t>::iterator start,
-                            std::vector<int32_t>& rpn)
-{
-  // parenthesis level counter
-  int counter = 0;
+std::vector<int32_t>::iterator
+CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
+                               const std::vector<int32_t>& rpn) {
 
-  // start at the beginning
+  int level = 0;
   auto it = start;
-  while(it < rpn.end()) {
+  while (it != rpn.begin()) {
+    // move through the rpn two at a time
     int32_t one = *it;
-    int32_t two = *(it+1);
+    int32_t two = *(it - 1);
 
+    // decrement parenthesis level if there are two adjacent surfaces
     if (one < OP_UNION && two < OP_UNION) {
-      counter++;
+      level--;
+    // increment if there are two adjacent operators
     } else if (one >= OP_UNION && two >= OP_UNION) {
-      counter--;
-      while(*it >= OP_UNION) { it++; }
-      it--;
+      level++;
     }
 
-    if (counter == 0) { break; }
-    it++;
+    // if the level gets to zero, return the position
+    if (level == 0) {
+      // move the iterator back one before leaving the loop
+      // so that all tokens in the parenthesis are included
+      it--;
+      break;
+    }
+
+    // continue loop
+    it--;
   }
   return it;
 }
 
-BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> 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
+void CSGCell::remove_complements(std::vector<int32_t>& rpn) {
   if (rpn.back() == OP_COMPLEMENT) {
-      rpn.pop_back();
-      apply_demorgan(rpn.begin(), rpn.end());
+    apply_demorgan(rpn.begin(), rpn.end());
+    rpn.pop_back();
   }
 
-  // check for parenthesis with complement at the start of the rpn
-  auto initial_right = find_right_parenthesis(rpn.begin(), rpn);
-  if (initial_right != rpn.end() && *initial_right == OP_COMPLEMENT) {
-    apply_demorgan(rpn.begin(), initial_right);
-    rpn.erase(initial_right);
+  std::vector<int32_t>::iterator 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);
+    // apply DeMorgan's law to any surfaces/operators between these
+    // positions in the RPN
+    apply_demorgan(left, it);
+    rpn.erase(it);
+
+    // update iterator position
+    it = std::find(rpn.begin(), rpn.end(), OP_COMPLEMENT);
   }
+}
+
+BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
+
+  remove_complements(rpn);
 
   auto it = rpn.begin();
   BoundingBox current = model::surfaces[abs(*it) - 1]->bounding_box(*it > 0);
@@ -659,10 +672,8 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
 
   while (it < rpn.end()) {
     // move through the rpn in twos
-    int32_t one = *it;
-    it++;
-    int32_t two = *it;
-    it++;
+    int32_t one = *it; it++;
+    int32_t two = *it; it++;
 
     // the first token should always be a surface
     Expects(one < OP_UNION);
@@ -682,19 +693,12 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
       // add until last two tokens in the sub-rpn are operators
       // (indicates a right parenthesis)
       while (!((*it >= OP_UNION) && (*(it + 1) >= OP_UNION))) {
-        subrpn.push_back(*it);
-        it++;
+        subrpn.push_back(*it); it++;
       }
 
       // add first operator to the subrpn
-      subrpn.push_back(*it);
-      it++;
+      subrpn.push_back(*it); it++;
 
-      // handle complement case using De Morgan's laws
-      if (*it == OP_COMPLEMENT) {
-        apply_demorgan(subrpn.begin(), subrpn.end());
-        it++;
-      }
       // save the last operator, tells us how to combine this region
       // with our current bounding box
       int32_t op = *it;

From c7471af7bfacbb742973ab5fa143678bcec3d888 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Wed, 11 Sep 2019 10:01:34 -0500
Subject: [PATCH 06/15] Clenaing and comments. No longer need to address the
 special case of a complement at the back of the RPN.

---
 src/cell.cpp | 14 +++++---------
 1 file changed, 5 insertions(+), 9 deletions(-)

diff --git a/src/cell.cpp b/src/cell.cpp
index f33926b303..69d398fda6 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -612,7 +612,7 @@ void CSGCell::apply_demorgan(std::vector<int32_t>::iterator start,
 std::vector<int32_t>::iterator
 CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
                                const std::vector<int32_t>& rpn) {
-
+  // start search at zero
   int level = 0;
   auto it = start;
   while (it != rpn.begin()) {
@@ -643,12 +643,7 @@ CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
 }
 
 void CSGCell::remove_complements(std::vector<int32_t>& rpn) {
-  if (rpn.back() == OP_COMPLEMENT) {
-    apply_demorgan(rpn.begin(), rpn.end());
-    rpn.pop_back();
-  }
-
-  std::vector<int32_t>::iterator it = std::find(rpn.begin(), rpn.end(), OP_COMPLEMENT);
+  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);
@@ -663,13 +658,15 @@ void CSGCell::remove_complements(std::vector<int32_t>& rpn) {
 }
 
 BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
-
+  // remove complements by adjusting surface signs and operators
   remove_complements(rpn);
 
+  // use the first token to set the bounding box
   auto it = rpn.begin();
   BoundingBox current = model::surfaces[abs(*it) - 1]->bounding_box(*it > 0);
   it++;
 
+  // loop over tokens
   while (it < rpn.end()) {
     // move through the rpn in twos
     int32_t one = *it; it++;
@@ -713,7 +710,6 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
       }
     }
   }
-
   return current;
 }
 

From 3dfc2da2a5e82ec12cc3cca0cb6f8493c7de71e2 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Wed, 11 Sep 2019 11:00:08 -0500
Subject: [PATCH 07/15] Removing unused bounding box method.

---
 include/openmc/surface.h | 7 -------
 1 file changed, 7 deletions(-)

diff --git a/include/openmc/surface.h b/include/openmc/surface.h
index 9cdc460bb9..d61ae07cef 100644
--- a/include/openmc/surface.h
+++ b/include/openmc/surface.h
@@ -83,13 +83,6 @@ struct BoundingBox
     return *this;
   }
 
-  inline bool is_valid() {
-    return xmin <= xmax &&
-           ymin <= ymax &&
-           zmin <= zmax;
-  }
-
-
 };
 
 //==============================================================================

From 5bb3d582c9c1ae815c332861ea5df0063b3ca337 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Thu, 12 Sep 2019 09:24:26 -0500
Subject: [PATCH 08/15] Consolidating increment operators. Using iterator to
 setup the subrpn to avoid eating memory if something goes wrong.

---
 src/cell.cpp | 26 +++++++++++++-------------
 1 file changed, 13 insertions(+), 13 deletions(-)

diff --git a/src/cell.cpp b/src/cell.cpp
index 69d398fda6..f560698aa2 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -616,7 +616,7 @@ CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
   int level = 0;
   auto it = start;
   while (it != rpn.begin()) {
-    // move through the rpn two at a time
+    // look at two tokens at a time
     int32_t one = *it;
     int32_t two = *(it - 1);
 
@@ -636,7 +636,7 @@ CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
       break;
     }
 
-    // continue loop
+    // continue loop, one token at a time
     it--;
   }
   return it;
@@ -669,8 +669,8 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
   // loop over tokens
   while (it < rpn.end()) {
     // move through the rpn in twos
-    int32_t one = *it; it++;
-    int32_t two = *it; it++;
+    int32_t one = *it++;
+    int32_t two = *it++;
 
     // the first token should always be a surface
     Expects(one < OP_UNION);
@@ -684,22 +684,22 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
     } else {
       // two surfaces in a row (left parenthesis),
       // create sub-rpn for region in parenthesis
-      std::vector<int32_t> 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 (!((*it >= OP_UNION) && (*(it + 1) >= OP_UNION))) {
-        subrpn.push_back(*it); it++;
+      auto subrpn_start = it - 2; // include current tokens
+      auto subrpn_end = it;
+      while (!((*subrpn_end >= OP_UNION) && (*(subrpn_end + 1) >= OP_UNION))) {
+        subrpn_end++;
       }
 
-      // add first operator to the subrpn
-      subrpn.push_back(*it); it++;
+      // create the subrpn, including the first of our current tokens
+      std::vector<int32_t> subrpn(subrpn_start, ++subrpn_end);
 
       // save the last operator, tells us how to combine this region
       // with our current bounding box
-      int32_t op = *it;
-      it++;
+      int32_t op = *subrpn_end++;
+      it = subrpn_end;
+
       // get bounding box for the subrpn
       BoundingBox sub_box = bounding_box_complex(subrpn);
       // combine the sub-rpn bounding box with our current cell box

From 25df513a77819d275cc5d98d2a1d5100f3c0c037 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Thu, 12 Sep 2019 09:29:12 -0500
Subject: [PATCH 09/15] Renaming function.

---
 include/openmc/cell.h | 2 +-
 src/cell.cpp          | 4 ++--
 2 files changed, 3 insertions(+), 3 deletions(-)

diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 90c03e3a90..2cd98ada1a 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -214,7 +214,7 @@ protected:
   static BoundingBox bounding_box_complex(std::vector<int32_t> rpn);
   static void apply_demorgan(std::vector<int32_t>::iterator start,
                              std::vector<int32_t>::iterator stop);
-  static void remove_complements(std::vector<int32_t>& rpn);
+  static void remove_complement_ops(std::vector<int32_t>& rpn);
   static std::vector<int32_t>::iterator
   find_left_parenthesis(std::vector<int32_t>::iterator start,
                         const std::vector<int32_t>& rpn);
diff --git a/src/cell.cpp b/src/cell.cpp
index f560698aa2..7ad48bd326 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -642,7 +642,7 @@ CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
   return it;
 }
 
-void CSGCell::remove_complements(std::vector<int32_t>& rpn) {
+void CSGCell::remove_complement_ops(std::vector<int32_t>& rpn) {
   auto it = std::find(rpn.begin(), rpn.end(), OP_COMPLEMENT);
   while (it != rpn.end()) {
     // find the opening parenthesis (if any)
@@ -659,7 +659,7 @@ void CSGCell::remove_complements(std::vector<int32_t>& rpn) {
 
 BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
   // remove complements by adjusting surface signs and operators
-  remove_complements(rpn);
+  remove_complement_ops(rpn);
 
   // use the first token to set the bounding box
   auto it = rpn.begin();

From 36fb89dead8d79b5b144d3035ea940ae6de9a719 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Fri, 13 Sep 2019 07:58:27 -0500
Subject: [PATCH 10/15] Making analagous function for finding the right
 parenthesis.

---
 include/openmc/cell.h |  4 ++++
 src/cell.cpp          | 43 +++++++++++++++++++++++++++++++++++++------
 2 files changed, 41 insertions(+), 6 deletions(-)

diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 2cd98ada1a..f6524c6ee3 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -218,6 +218,10 @@ protected:
   static std::vector<int32_t>::iterator
   find_left_parenthesis(std::vector<int32_t>::iterator start,
                         const std::vector<int32_t>& rpn);
+  static std::vector<int32_t>::iterator
+  find_right_parenthesis(std::vector<int32_t>::iterator start,
+                        const std::vector<int32_t>& rpn);
+
 };
 
 //==============================================================================
diff --git a/src/cell.cpp b/src/cell.cpp
index 7ad48bd326..bc0acbe73d 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -642,6 +642,40 @@ CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
   return it;
 }
 
+std::vector<int32_t>::iterator
+CSGCell::find_right_parenthesis(std::vector<int32_t>::iterator start,
+                                const std::vector<int32_t>& rpn) {
+  // start search at zero
+  int level = 0;
+  auto it = start;
+  while (it != rpn.end()) {
+    // 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) {
+      level++;
+    // increment if there are two adjacent operators
+    } else if (one >= OP_UNION && two >= OP_UNION) {
+      level--;
+    }
+
+    // if the level gets to zero, return the position
+    if (level == 0) {
+      // move the iterator back one before leaving the loop
+      // so that all tokens in the parenthesis are included
+      it++;
+      break;
+    }
+
+    // continue loop, one token at a time
+    it++;
+  }
+  return it;
+}
+
+
 void CSGCell::remove_complement_ops(std::vector<int32_t>& rpn) {
   auto it = std::find(rpn.begin(), rpn.end(), OP_COMPLEMENT);
   while (it != rpn.end()) {
@@ -687,20 +721,17 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
       // add until last two tokens in the sub-rpn are operators
       // (indicates a right parenthesis)
       auto subrpn_start = it - 2; // include current tokens
-      auto subrpn_end = it;
-      while (!((*subrpn_end >= OP_UNION) && (*(subrpn_end + 1) >= OP_UNION))) {
-        subrpn_end++;
-      }
+      auto subrpn_end = find_right_parenthesis(subrpn_start, rpn);
 
       // create the subrpn, including the first of our current tokens
-      std::vector<int32_t> subrpn(subrpn_start, ++subrpn_end);
+      std::vector<int32_t> subrpn(subrpn_start, subrpn_end);
 
       // save the last operator, tells us how to combine this region
       // with our current bounding box
       int32_t op = *subrpn_end++;
       it = subrpn_end;
 
-      // get bounding box for the subrpn
+      // get bounding box for the sub-rpn
       BoundingBox sub_box = bounding_box_complex(subrpn);
       // combine the sub-rpn bounding box with our current cell box
       if (op == OP_UNION) {

From 962f3d46f1947d3d8584f66dd9770c9a41bc1259 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Fri, 13 Sep 2019 08:29:22 -0500
Subject: [PATCH 11/15] Adding documentation for some of these new methods.

---
 include/openmc/cell.h | 12 ++++++++++++
 1 file changed, 12 insertions(+)

diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index f6524c6ee3..5e53e3bc43 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -212,12 +212,24 @@ 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<int32_t> rpn);
+
+  //! Applies DeMorgan's laws to a section of the RPN
   static void apply_demorgan(std::vector<int32_t>::iterator start,
                              std::vector<int32_t>::iterator stop);
+
+  //! Removes complement operators from the RPN
   static void remove_complement_ops(std::vector<int32_t>& 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)
   static std::vector<int32_t>::iterator
   find_left_parenthesis(std::vector<int32_t>::iterator start,
                         const std::vector<int32_t>& rpn);
+
+  //! Returns the ending position of a parenthesis block (immediately after two
+  //! operator tokens) in the RPN given a starting position at the beginning of
+  //! that block (immediately before two surface tokens)
   static std::vector<int32_t>::iterator
   find_right_parenthesis(std::vector<int32_t>::iterator start,
                         const std::vector<int32_t>& rpn);

From ec266e3e0b777bbe8a8e0238bf3809f8e34a60b7 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Fri, 13 Sep 2019 11:35:20 -0500
Subject: [PATCH 12/15] Updating comments and a variable name

---
 src/cell.cpp | 34 ++++++++++++++++------------------
 1 file changed, 16 insertions(+), 18 deletions(-)

diff --git a/src/cell.cpp b/src/cell.cpp
index bc0acbe73d..7beb98cc2f 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -613,7 +613,7 @@ std::vector<int32_t>::iterator
 CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
                                const std::vector<int32_t>& rpn) {
   // start search at zero
-  int level = 0;
+  int parenthesis_level = 0;
   auto it = start;
   while (it != rpn.begin()) {
     // look at two tokens at a time
@@ -622,16 +622,16 @@ CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
 
     // decrement parenthesis level if there are two adjacent surfaces
     if (one < OP_UNION && two < OP_UNION) {
-      level--;
+      parenthesis_level--;
     // increment if there are two adjacent operators
     } else if (one >= OP_UNION && two >= OP_UNION) {
-      level++;
+      parenthesis_level++;
     }
 
     // if the level gets to zero, return the position
-    if (level == 0) {
+    if (parenthesis_level == 0) {
       // move the iterator back one before leaving the loop
-      // so that all tokens in the parenthesis are included
+      // so that all tokens in the parenthesis block are included
       it--;
       break;
     }
@@ -645,26 +645,25 @@ CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
 std::vector<int32_t>::iterator
 CSGCell::find_right_parenthesis(std::vector<int32_t>::iterator start,
                                 const std::vector<int32_t>& rpn) {
-  // start search at zero
-  int level = 0;
+  int parenthesis_level = 0;
   auto it = start;
   while (it != rpn.end()) {
     // 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
+    // increment parenthesis level if there are two adjacent surfaces
     if (one < OP_UNION && two < OP_UNION) {
-      level++;
-    // increment if there are two adjacent operators
+      parenthesis_level++;
+    // decrement if there are two adjacent operators
     } else if (one >= OP_UNION && two >= OP_UNION) {
-      level--;
+      parenthesis_level--;
     }
 
     // if the level gets to zero, return the position
-    if (level == 0) {
-      // move the iterator back one before leaving the loop
-      // so that all tokens in the parenthesis are included
+    if (parenthesis_level == 0) {
+      // move the iterator forward one before leaving the loop
+      // so that all tokens in the parenthesis block are included
       it++;
       break;
     }
@@ -675,17 +674,18 @@ CSGCell::find_right_parenthesis(std::vector<int32_t>::iterator start,
   return it;
 }
 
-
 void CSGCell::remove_complement_ops(std::vector<int32_t>& 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<int32_t> 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);
   }
@@ -718,8 +718,6 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
     } else {
       // two surfaces in a row (left parenthesis),
       // create sub-rpn for region in parenthesis
-      // add until last two tokens in the sub-rpn are operators
-      // (indicates a right parenthesis)
       auto subrpn_start = it - 2; // include current tokens
       auto subrpn_end = find_right_parenthesis(subrpn_start, rpn);
 

From 8ebfdc3eb78ce68f7441c70f2c649ffe501c8e0e Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Fri, 13 Sep 2019 13:45:18 -0500
Subject: [PATCH 13/15] Improving documentation.

---
 include/openmc/cell.h | 7 +++++++
 src/cell.cpp          | 2 +-
 2 files changed, 8 insertions(+), 1 deletion(-)

diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 5e53e3bc43..5ba714570a 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -214,15 +214,20 @@ protected:
   static BoundingBox bounding_box_complex(std::vector<int32_t> 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<int32_t>::iterator start,
                              std::vector<int32_t>::iterator stop);
 
   //! Removes complement operators from the RPN
+  //! \param rpn The rpn to remove complement operators from.
   static void remove_complement_ops(std::vector<int32_t>& 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<int32_t>::iterator
   find_left_parenthesis(std::vector<int32_t>::iterator start,
                         const std::vector<int32_t>& rpn);
@@ -230,6 +235,8 @@ protected:
   //! Returns the ending position of a parenthesis block (immediately after two
   //! operator tokens) in the RPN given a starting position at the beginning of
   //! that block (immediately before two surface tokens)
+  //! \param start Starting position of the search
+  //! \param rpn The rpn being searched
   static std::vector<int32_t>::iterator
   find_right_parenthesis(std::vector<int32_t>::iterator start,
                         const std::vector<int32_t>& rpn);
diff --git a/src/cell.cpp b/src/cell.cpp
index 7beb98cc2f..d87af54e47 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -601,7 +601,7 @@ BoundingBox CSGCell::bounding_box_simple() const {
 void CSGCell::apply_demorgan(std::vector<int32_t>::iterator start,
                              std::vector<int32_t>::iterator stop)
 {
-  while(start < 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; }

From db2cf9eea693eb8725cff25db9fccafd6e5899b0 Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Mon, 16 Sep 2019 11:25:59 -0500
Subject: [PATCH 14/15] Updating bounding box to a stack-based traversal of the
 RPN for complex cells.

---
 src/cell.cpp | 59 ++++++++++++++++------------------------------------
 1 file changed, 18 insertions(+), 41 deletions(-)

diff --git a/src/cell.cpp b/src/cell.cpp
index d87af54e47..914391ae22 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -695,51 +695,28 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
   // remove complements by adjusting surface signs and operators
   remove_complement_ops(rpn);
 
-  // use the first token to set the bounding box
-  auto it = rpn.begin();
-  BoundingBox current = model::surfaces[abs(*it) - 1]->bounding_box(*it > 0);
-  it++;
+  std::vector<BoundingBox> stack(rpn.size());
+  int i_stack = -1;
 
-  // loop over tokens
-  while (it < rpn.end()) {
-    // move through the rpn in twos
-    int32_t one = *it++;
-    int32_t two = *it++;
-
-    // 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
-      auto subrpn_start = it - 2; // include current tokens
-      auto subrpn_end = find_right_parenthesis(subrpn_start, rpn);
-
-      // create the subrpn, including the first of our current tokens
-      std::vector<int32_t> subrpn(subrpn_start, subrpn_end);
-
-      // save the last operator, tells us how to combine this region
-      // with our current bounding box
-      int32_t op = *subrpn_end++;
-      it = subrpn_end;
-
-      // get bounding box for the sub-rpn
-      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;
+
+  if (i_stack == 0) {
+    return stack[i_stack];
+  } else {
+    return BoundingBox();
+  }
+
 }
 
 BoundingBox CSGCell::bounding_box() const {

From c0e32b751bd6b3d3e64939a43227adb7de079b4e Mon Sep 17 00:00:00 2001
From: Patrick Shriwise <pshriwise@gmail.com>
Date: Mon, 16 Sep 2019 23:48:52 -0500
Subject: [PATCH 15/15] Addressing comments from @paulromano.

---
 include/openmc/cell.h |  9 ---------
 src/cell.cpp          | 40 ++--------------------------------------
 2 files changed, 2 insertions(+), 47 deletions(-)

diff --git a/include/openmc/cell.h b/include/openmc/cell.h
index 5ba714570a..7206938c40 100644
--- a/include/openmc/cell.h
+++ b/include/openmc/cell.h
@@ -232,15 +232,6 @@ protected:
   find_left_parenthesis(std::vector<int32_t>::iterator start,
                         const std::vector<int32_t>& rpn);
 
-  //! Returns the ending position of a parenthesis block (immediately after two
-  //! operator tokens) in the RPN given a starting position at the beginning of
-  //! that block (immediately before two surface tokens)
-  //! \param start Starting position of the search
-  //! \param rpn The rpn being searched
-  static std::vector<int32_t>::iterator
-  find_right_parenthesis(std::vector<int32_t>::iterator start,
-                        const std::vector<int32_t>& rpn);
-
 };
 
 //==============================================================================
diff --git a/src/cell.cpp b/src/cell.cpp
index 914391ae22..3a7a2598fc 100644
--- a/src/cell.cpp
+++ b/src/cell.cpp
@@ -642,38 +642,6 @@ CSGCell::find_left_parenthesis(std::vector<int32_t>::iterator start,
   return it;
 }
 
-std::vector<int32_t>::iterator
-CSGCell::find_right_parenthesis(std::vector<int32_t>::iterator start,
-                                const std::vector<int32_t>& rpn) {
-  int parenthesis_level = 0;
-  auto it = start;
-  while (it != rpn.end()) {
-    // look at two tokens at a time
-    int32_t one = *it;
-    int32_t two = *(it + 1);
-
-    // increment parenthesis level if there are two adjacent surfaces
-    if (one < OP_UNION && two < OP_UNION) {
-      parenthesis_level++;
-    // decrement 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 forward 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<int32_t>& rpn) {
   auto it = std::find(rpn.begin(), rpn.end(), OP_COMPLEMENT);
   while (it != rpn.end()) {
@@ -711,12 +679,8 @@ BoundingBox CSGCell::bounding_box_complex(std::vector<int32_t> rpn) {
     }
   }
 
-  if (i_stack == 0) {
-    return stack[i_stack];
-  } else {
-    return BoundingBox();
-  }
-
+  Ensures(i_stack == 0);
+  return stack.front();
 }
 
 BoundingBox CSGCell::bounding_box() const {