Represent cell regions as expression trees

include/openmc/cell.h

@@ -66,12 +66,12 @@ class Region {
   //! \brief Determine if a cell contains the particle at a given location.
   //!
   //! The bounds of the cell are determined by a logical expression involving
-  //! surface half-spaces. The expression used is given in infix notation
+  //! surface half-spaces, stored as an expression tree of intersections and
+  //! unions with half-spaces as leaves.
   //!
   //! The function is split into two cases, one for simple cells (those
   //! involving only the intersection of half-spaces) and one for complex cells.
-  //! Both cases use short circuiting; however, in the case fo complex cells,
-  //! the complexity increases with the binary operators involved.
+  //! Both cases use short circuiting.
   //! \param r The 3D Cartesian coordinate to check.
   //! \param u A direction used to "break ties" the coordinates are very
   //!   close to a surface.
@@ -85,39 +85,50 @@ class Region {
     Position r, Direction u, int32_t on_surface) const;
 
   //! Get the BoundingBox for this cell.
-  BoundingBox bounding_box(int32_t cell_id) const;
+  BoundingBox bounding_box() const;
 
   //! Get the CSG expression as a string
   std::string str() const;
 
-  //! Get a vector containing all the surfaces in the region expression
+  //! Get a vector containing all the half-spaces in the region expression
   vector<int32_t> surfaces() const;
 
-  //! Get size of surfaces
-  int n_surfaces() const { return expression_.size(); }
+  //! Get the number of half-spaces in the region expression
+  int n_surfaces() const;
 
   //----------------------------------------------------------------------------
   // Accessors
 
   //! Get Boolean of if the cell is simple or not
-  bool is_simple() const { return simple_; }
+  bool is_simple() const { return !complex_; }
 
 private:
   //----------------------------------------------------------------------------
-  // Private Methods
+  // Types
+
+  //! Node of the region expression tree. Nodes are stored in pre-order, so
+  //! the children of an operator node follow it, and the subtree of a node
+  //! ends just before index end. Children of an operator node are never
+  //! operator nodes of the same type.
+  struct Node {
+    enum class Type : int8_t { HALFSPACE, INTERSECTION, UNION };
+    Type type;
+    int32_t halfspace; //!< Signed surface index + 1 for HALFSPACE nodes
+    int32_t end;       //!< Index one past the last node of the subtree
+    int32_t parent;    //!< Index of the parent node (-1 for the root)
+  };
 
-  //! Get a vector of the region expression in postfix notation
-  vector<int32_t> generate_postfix(int32_t cell_id) const;
+  //----------------------------------------------------------------------------
+  // Private Methods
 
   //! Determine if a particle is inside the cell for a simple cell (only
   //! intersection operators)
   bool contains_simple(Position r, Direction u, int32_t on_surface) const;
 
   //! Determine if a particle is inside the cell for a complex cell.
   //!
-  //! Uses the combination of half-spaces and binary operators to determine
-  //! if short circuiting can be used. Short circuiting uses the relative and
-  //! absolute depth of parentheses in the expression.
+  //! Evaluates the expression tree, skipping the remaining children of an
+  //! operator node as soon as its value is known.
   bool contains_complex(Position r, Direction u, int32_t on_surface) const;
 
   //! Find the nearest intersection with any surface in the region expression.
@@ -128,32 +139,21 @@ class Region {
   std::pair<double, int32_t> distance_complex(
     Position r, Direction u, int32_t on_surface) const;
 
-  //! BoundingBox if the particle is in a simple cell.
-  BoundingBox bounding_box_simple() const;
-
-  //! BoundingBox if the particle is in a complex cell.
-  BoundingBox bounding_box_complex(vector<int32_t> postfix) const;
-
-  //! Enforce precedence between intersections and unions
-  void enforce_precedence();
-
-  //! Add parenthesis to enforce precedence
-  void add_parentheses(int64_t start);
-
-  //! Remove complement operators from the expression
-  void remove_complement_ops();
-
-  //! Remove complement operators by using DeMorgan's laws
-  void apply_demorgan(
-    vector<int32_t>::iterator start, vector<int32_t>::iterator stop);
-
   //----------------------------------------------------------------------------
   // Private Data
 
-  //! Definition of spatial region as Boolean expression of half-spaces
-  // TODO: Should this be a vector of some other type
-  vector<int32_t> expression_;
-  bool simple_; //!< Does the region contain only intersections?
+  //! Signed surface indices + 1 of the half-spaces in the region expression,
+  //! in order. A simple region is the intersection of these half-spaces.
+  vector<int32_t> halfspaces_;
+
+  //! Data needed only by complex regions, kept out of line so that regions,
+  //! and the cells holding them, stay small for simple cells
+  struct Complex {
+    vector<Node> nodes; //!< Expression tree in pre-order
+  };
+
+  //! Data of a complex region (null for a simple region)
+  unique_ptr<Complex> complex_;
 };
 
 //==============================================================================
@@ -447,26 +447,14 @@ class CSGCell : public Cell {
     return region_.contains(r, u, on_surface);
   }
 
-  BoundingBox bounding_box() const override
-  {
-    return region_.bounding_box(id_);
-  }
+  BoundingBox bounding_box() const override { return region_.bounding_box(); }
 
   void to_hdf5_inner(hid_t group_id) const override;
 
   bool is_simple() const override { return region_.is_simple(); }
 
   virtual GeometryType geom_type() const override { return GeometryType::CSG; }
 
-protected:
-  //! Returns the beginning position of a parenthesis block (immediately before
-  //! two surface tokens) in the RPN given a starting position at the end of
-  //! that block (immediately after two surface tokens)
-  //! \param start Starting position of the search
-  //! \param rpn The rpn being searched
-  static vector<int32_t>::iterator find_left_parenthesis(
-    vector<int32_t>::iterator start, const vector<int32_t>& rpn);
-
 private:
   Region region_;
 };