Source code for openmc.region

from abc import ABC, abstractmethod
from collections import OrderedDict
from collections.abc import MutableSequence
from copy import deepcopy

import numpy as np

from .checkvalue import check_type


[docs]class Region(ABC): """Region of space that can be assigned to a cell. Region is an abstract base class that is inherited by :class:`openmc.Halfspace`, :class:`openmc.Intersection`, :class:`openmc.Union`, and :class:`openmc.Complement`. Each of those respective classes are typically not instantiated directly but rather are created through operators of the Surface and Region classes. """ def __and__(self, other): return Intersection((self, other)) def __or__(self, other): return Union((self, other)) def __invert__(self): return Complement(self) @abstractmethod def __contains__(self, point): pass @abstractmethod def __str__(self): pass def __eq__(self, other): if not isinstance(other, type(self)): return False else: return str(self) == str(other)
[docs] def get_surfaces(self, surfaces=None): """Recursively find all surfaces referenced by a region and return them Parameters ---------- surfaces: collections.OrderedDict, optional Dictionary mapping surface IDs to :class:`openmc.Surface` instances Returns ------- surfaces: collections.OrderedDict Dictionary mapping surface IDs to :class:`openmc.Surface` instances """ if surfaces is None: surfaces = OrderedDict() for region in self: surfaces = region.get_surfaces(surfaces) return surfaces
[docs] def remove_redundant_surfaces(self, redundant_surfaces): """Recursively remove all redundant surfaces referenced by this region .. versionadded:: 0.12 Parameters ---------- redundant_surfaces : dict Dictionary mapping redundant surface IDs to class:`openmc.Surface` instances that should replace them. """ for region in self: region.remove_redundant_surfaces(redundant_surfaces)
[docs] @staticmethod def from_expression(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 surface IDs that appear in the Boolean expression and whose values are Surface objects. """ # Strip leading and trailing whitespace expression = expression.strip() # 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 appropriate half-space if i_start >= 0: # When an opening parenthesis appears after a non-operator, # there's an implicit intersection operator between them if expression[i] == '(': tokens.append(' ') j = int(expression[i_start:i]) if j < 0: tokens.append(-surfaces[abs(j)]) else: tokens.append(+surfaces[abs(j)]) if expression[i] in '()|~': # For everything other than intersection, add the operator # to the list of tokens tokens.append(expression[i]) # If two parentheses appear immediately adjacent to one # another, we need an intersection between them if expression[i:i+2] == ')(': tokens.append(' ') else: # Find next non-space character while expression[i+1] == ' ': i += 1 # If previous token is a halfspace or right parenthesis and # next token is not a left parenthesis or union operator, # that implies that the whitespace is to be interpreted as # an intersection operator if (i_start >= 0 or tokens[-1] == ')') and \ expression[i+1] not in ')|': tokens.append(' ') i_start = -1 else: # Check for invalid characters if expression[i] not in '-+0123456789': raise SyntaxError(f"Invalid character '{expression[i]}' 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)]) else: tokens.append(+surfaces[abs(j)]) # The functions below are used to apply an operator to operands on the # output queue during the shunting yard algorithm. def can_be_combined(region): return isinstance(region, Complement) or hasattr(region, 'surface') def apply_operator(output, operator): r2 = output.pop() if operator == ' ': r1 = output.pop() if isinstance(r1, Intersection): r1 &= r2 output.append(r1) elif isinstance(r2, Intersection) and can_be_combined(r1): r2.insert(0, r1) output.append(r2) else: output.append(r1 & r2) elif operator == '|': r1 = output.pop() if isinstance(r1, Union): r1 |= r2 output.append(r1) elif isinstance(r2, Union) and can_be_combined(r1): r2.insert(0, r1) output.append(r2) else: output.append(r1 | r2) elif operator == '~': output.append(~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]
[docs] def clone(self, memo=None): """Create a copy of this region - each of the surfaces in the region's nodes will be cloned and will have new unique IDs. Parameters ---------- memo : dict or None A nested dictionary of previously cloned objects. This parameter is used internally and should not be specified by the user. Returns ------- clone : openmc.Region The clone of this region """ if memo is None: memo = {} clone = deepcopy(self) clone[:] = [n.clone(memo) for n in self] return clone
[docs] def translate(self, vector, inplace=False, memo=None): """Translate region in given direction Parameters ---------- vector : iterable of float Direction in which region should be translated inplace : bool Whether or not to return a region based on new surfaces or one based on the original surfaces that have been modified. memo : dict or None Dictionary used for memoization. This parameter is used internally and should not be specified by the user. Returns ------- openmc.Region Translated region """ if memo is None: memo = {} return type(self)(n.translate(vector, inplace, memo) for n in self)
[docs] def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False, memo=None): r"""Rotate surface by angles provided or by applying matrix directly. .. versionadded:: 0.12 Parameters ---------- rotation : 3-tuple of float, or 3x3 iterable A 3-tuple of angles :math:`(\phi, \theta, \psi)` in degrees where the first element is the rotation about the x-axis in the fixed laboratory frame, the second element is the rotation about the y-axis in the fixed laboratory frame, and the third element is the rotation about the z-axis in the fixed laboratory frame. The rotations are active rotations. Additionally a 3x3 rotation matrix can be specified directly either as a nested iterable or array. pivot : iterable of float, optional (x, y, z) coordinates for the point to rotate about. Defaults to (0., 0., 0.) order : str, optional A string of 'x', 'y', and 'z' in some order specifying which rotation to perform first, second, and third. Defaults to 'xyz' which means, the rotation by angle :math:`\phi` about x will be applied first, followed by :math:`\theta` about y and then :math:`\psi` about z. This corresponds to an x-y-z extrinsic rotation as well as a z-y'-x'' intrinsic rotation using Tait-Bryan angles :math:`(\phi, \theta, \psi)`. inplace : bool Whether or not to return a new instance of Surface or to modify the coefficients of this Surface in place. Defaults to False. memo : dict or None Dictionary used for memoization Returns ------- openmc.Region Translated region """ if memo is None: memo = {} return type(self)(n.rotate(rotation, pivot=pivot, order=order, inplace=inplace, memo=memo) for n in self)
[docs]class Intersection(Region, MutableSequence): r"""Intersection of two or more regions. Instances of Intersection are generally created via the & operator applied to two instances of :class:`openmc.Region`. This is illustrated in the following example: >>> equator = openmc.ZPlane(z0=0.0) >>> earth = openmc.Sphere(r=637.1e6) >>> northern_hemisphere = -earth & +equator >>> southern_hemisphere = -earth & -equator >>> type(northern_hemisphere) <class 'openmc.region.Intersection'> Instances of this class behave like a mutable sequence, e.g., they can be indexed and have an append() method. Parameters ---------- nodes : iterable of openmc.Region Regions to take the intersection of Attributes ---------- bounding_box : tuple of numpy.array Lower-left and upper-right coordinates of an axis-aligned bounding box """ def __init__(self, nodes): self._nodes = list(nodes) def __and__(self, other): new = Intersection(self) new &= other return new def __iand__(self, other): if isinstance(other, Intersection): self.extend(other) else: self.append(other) return self # Implement mutable sequence protocol by delegating to list def __getitem__(self, key): return self._nodes[key] def __setitem__(self, key, value): self._nodes[key] = value def __delitem__(self, key): del self._nodes[key] def __len__(self): return len(self._nodes)
[docs] def insert(self, index, value): self._nodes.insert(index, value)
def __contains__(self, point): """Check whether a point is contained in the region. Parameters ---------- point : 3-tuple of float Cartesian coordinates, :math:`(x',y',z')`, of the point Returns ------- bool Whether the point is in the region """ return all(point in n for n in self) def __str__(self): return '(' + ' '.join(map(str, self)) + ')' @property def bounding_box(self): lower_left = np.array([-np.inf, -np.inf, -np.inf]) upper_right = np.array([np.inf, np.inf, np.inf]) for n in self: lower_left_n, upper_right_n = n.bounding_box lower_left[:] = np.maximum(lower_left, lower_left_n) upper_right[:] = np.minimum(upper_right, upper_right_n) return lower_left, upper_right
[docs]class Union(Region, MutableSequence): r"""Union of two or more regions. Instances of Union are generally created via the | operator applied to two instances of :class:`openmc.Region`. This is illustrated in the following example: >>> s1 = openmc.ZPlane(z0=0.0) >>> s2 = openmc.Sphere(r=637.1e6) >>> type(-s2 | +s1) <class 'openmc.region.Union'> Instances of this class behave like a mutable sequence, e.g., they can be indexed and have an append() method. Parameters ---------- nodes : iterable of openmc.Region Regions to take the union of Attributes ---------- bounding_box : 2-tuple of numpy.array Lower-left and upper-right coordinates of an axis-aligned bounding box """ def __init__(self, nodes): self._nodes = list(nodes) def __or__(self, other): new = Union(self) new |= other return new def __ior__(self, other): if isinstance(other, Union): self.extend(other) else: self.append(other) return self # Implement mutable sequence protocol by delegating to list def __getitem__(self, key): return self._nodes[key] def __setitem__(self, key, value): self._nodes[key] = value def __delitem__(self, key): del self._nodes[key] def __len__(self): return len(self._nodes)
[docs] def insert(self, index, value): self._nodes.insert(index, value)
def __contains__(self, point): """Check whether a point is contained in the region. Parameters ---------- point : 3-tuple of float Cartesian coordinates, :math:`(x',y',z')`, of the point Returns ------- bool Whether the point is in the region """ return any(point in n for n in self) def __str__(self): return '(' + ' | '.join(map(str, self)) + ')' @property def bounding_box(self): lower_left = np.array([np.inf, np.inf, np.inf]) upper_right = np.array([-np.inf, -np.inf, -np.inf]) for n in self: lower_left_n, upper_right_n = n.bounding_box lower_left[:] = np.minimum(lower_left, lower_left_n) upper_right[:] = np.maximum(upper_right, upper_right_n) return lower_left, upper_right
[docs]class Complement(Region): """Complement of a region. The Complement of an existing :class:`openmc.Region` can be created by using the ~ operator as the following example demonstrates: >>> xl = openmc.XPlane(-10.0) >>> xr = openmc.XPlane(10.0) >>> yl = openmc.YPlane(-10.0) >>> yr = openmc.YPlane(10.0) >>> inside_box = +xl & -xr & +yl & -yr >>> outside_box = ~inside_box >>> type(outside_box) <class 'openmc.region.Complement'> Parameters ---------- node : openmc.Region Region to take the complement of Attributes ---------- node : openmc.Region Regions to take the complement of bounding_box : tuple of numpy.array Lower-left and upper-right coordinates of an axis-aligned bounding box """ def __init__(self, node): self.node = node def __contains__(self, point): """Check whether a point is contained in the region. Parameters ---------- point : 3-tuple of float Cartesian coordinates, :math:`(x',y',z')`, of the point Returns ------- bool Whether the point is in the region """ return point not in self.node def __str__(self): return '~' + str(self.node) @property def node(self): return self._node @node.setter def node(self, node): check_type('node', node, Region) self._node = node @property def bounding_box(self): # Use De Morgan's laws to distribute the complement operator so that it # only applies to surface half-spaces, thus allowing us to calculate the # bounding box in the usual recursive manner. if isinstance(self.node, Union): temp_region = Intersection(~n for n in self.node) elif isinstance(self.node, Intersection): temp_region = Union(~n for n in self.node) elif isinstance(self.node, Complement): temp_region = self.node.node else: temp_region = ~self.node return temp_region.bounding_box
[docs] def get_surfaces(self, surfaces=None): """Recursively find and return all the surfaces referenced by the node Parameters ---------- surfaces: collections.OrderedDict, optional Dictionary mapping surface IDs to :class:`openmc.Surface` instances Returns ------- surfaces: collections.OrderedDict Dictionary mapping surface IDs to :class:`openmc.Surface` instances """ if surfaces is None: surfaces = OrderedDict() for region in self.node: surfaces = region.get_surfaces(surfaces) return surfaces
[docs] def remove_redundant_surfaces(self, redundant_surfaces): """Recursively remove all redundant surfaces referenced by this region .. versionadded:: 0.12 Parameters ---------- redundant_surfaces : dict Dictionary mapping redundant surface IDs to class:`openmc.Surface` instances that should replace them. """ for region in self.node: region.remove_redundant_surfaces(redundant_surfaces)
[docs] def clone(self, memo=None): if memo is None: memo = {} clone = deepcopy(self) clone.node = self.node.clone(memo) return clone
[docs] def translate(self, vector, inplace=False, memo=None): if memo is None: memo = {} return type(self)(self.node.translate(vector, inplace, memo))
[docs] def rotate(self, rotation, pivot=(0., 0., 0.), order='xyz', inplace=False, memo=None): if memo is None: memo = {} return type(self)(self.node.rotate(rotation, pivot=pivot, order=order, inplace=inplace, memo=memo))