Add Region.from_expression method

openmc/region.py

@@ -11,6 +11,132 @@ class Region(object):
     def __str__(self):
         return ''
 
+    @classmethod
+    def from_expression(cls, expression, surfaces):
+        """Generate a region given an infix expression.
+
+        Parameters
+        ----------
+        expression : str
+            Boolean expression relating surface half-spaces. The possible
+            operators are union '^', intersection ' ', and complement '~'. For
+            example, '(1 -2) ^ 3 ~(4 -5)'.
+        surfaces : dict
+            Dictionary whose keys are suface IDs that appear in the Boolean
+            expression and whose values are Surface objects.
+
+        """
+
+        # Convert the string expression into a list of tokens, i.e., operators
+        # and surface half-spaces, representing the expression in infix
+        # notation.
+        i = 0
+        i_start = -1
+        tokens = []
+        while i < len(expression):
+            if expression[i] in '()^~ ':
+                # If special character appears immediately after a non-operator,
+                # create a token with the apporpriate half-space
+                if i_start >= 0:
+                    j = int(expression[i_start:i])
+                    if j < 0:
+                        tokens.append(surfaces[abs(j)].negative)
+                    else:
+                        tokens.append(surfaces[abs(j)].positive)
+
+                if expression[i] in '()^~':
+                    # For everything other than intersection, add the operator
+                    # to the list of tokens
+                    tokens.append(expression[i])
+                else:
+                    # For spaces, we need to check the context further. If it
+                    # doesn't appear before a right parentheses or union, it is
+                    # interpreted to be as an intersection operator
+                    while expression[i+1] == ' ':
+                        i += 1
+
+                    if i_start >= 0 and expression[i+1] not in ')^':
+                        tokens.append(' ')
+
+                i_start = -1
+            else:
+                # Check for invalid characters
+                if expression[i] not in '-0123456789':
+                    raise SyntaxError('Invalid character in expression')
+
+                # If we haven't yet reached the start of a word, start one
+                if i_start < 0:
+                    i_start = i
+            i += 1
+
+        # If we've reached the end and we're still in a word, create a
+        # half-space token and add it to the list
+        if i_start >= 0:
+            j = int(expression[i_start:])
+            if j < 0:
+                tokens.append(surfaces[abs(j)].negative)
+            else:
+                tokens.append(surfaces[abs(j)].positive)
+
+        # This function is used below to apply an operator to operands on the
+        # output queue during the shunting yard algorithm.
+        def apply_operator(output, operator):
+            r2 = output.pop()
+            if operator == ' ':
+                r1 = output.pop()
+                output.append(Intersection(r1, r2))
+            elif operator == '^':
+                r1 = output.pop()
+                output.append(Union(r1, r2))
+            elif operator == '~':
+                output.append(Complement(r2))
+
+
+        # The following is an implementation of the shunting yard algorithm to
+        # generate an abstract syntax tree for the region expression.
+        output = []
+        stack = []
+        precedence = {'^': 1, ' ': 2, '~': 3}
+        associativity = {'^': 'left', ' ': 'left', '~': 'right'}
+        for token in tokens:
+            if token in (' ', '^', '~'):
+                # Normal operators
+                while stack:
+                    op = stack[-1]
+                    if (op not in ('(', ')') and
+                        ((associativity[token] == 'right' and
+                          precedence[token] < precedence[op]) or
+                         (associativity[token] == 'left' and
+                          precedence[token] <= precedence[op]))):
+                        apply_operator(output, stack.pop())
+                    else:
+                        break
+                stack.append(token)
+            elif token == '(':
+                # Left parentheses
+                stack.append(token)
+            elif token == ')':
+                # Right parentheses
+                while stack[-1] != '(':
+                    apply_operator(output, stack.pop())
+                    if len(stack) == 0:
+                        raise SyntaxError('Mismatched parentheses in '
+                                          'region specification.')
+                stack.pop()
+            else:
+                # Surface halfspaces
+                output.append(token)
+        while stack:
+            if stack[-1] in '()':
+                raise SyntaxError('Mismatched parentheses in region '
+                                  'specification.')
+            apply_operator(output, stack.pop())
+
+        # Since we are generating an abstract syntax tree rather than a reverse
+        # Polish notation expression, the output queue should have a single item
+        # at the end
+        return output[0]
+
 
 class Intersection(Region):
     """Intersection of two or more regions.

openmc/summary.py

@@ -1,6 +1,7 @@
 import numpy as np
 
 import openmc
+from openmc.region import Region
 
 
 class Summary(object):
@@ -240,16 +241,10 @@ class Summary(object):
             # Store Cell fill information for after Universe/Lattice creation
             self._cell_fills[index] = (fill_type, fill)
 
-            # Iterate over all Surfaces and add them to the Cell
-            for token in region.split():
-                try:
-                    surface_halfspace = int(token)
-                    halfspace = np.sign(surface_halfspace)
-                    surface_id = abs(surface_halfspace)
-                    surface = self.get_surface_by_id(surface_id)
-                    cell.add_surface(surface, halfspace)
-                except ValueError:
-                    pass
+            # Generate Region object given infix expression
+            if region:
+                cell.region = Region.from_expression(
+                    region, {s.id: s for s in self.surfaces.values()})
 
             # Add the Cell to the global dictionary of all Cells
             self.cells[index] = cell