diff --git a/CMakeLists.txt b/CMakeLists.txt index d628cd2a6e..5e74ce96dd 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -9,6 +9,11 @@ set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${CMAKE_BINARY_DIR}/bin) # Set module path set(CMAKE_MODULE_PATH ${CMAKE_CURRENT_SOURCE_DIR}/cmake/Modules) +# Allow user to specify _ROOT variables +if (NOT (CMAKE_VERSION VERSION_LESS 3.12)) + cmake_policy(SET CMP0074 NEW) +endif() + #=============================================================================== # Command line options #=============================================================================== @@ -35,18 +40,18 @@ endif() #=============================================================================== if(dagmc) find_package(DAGMC REQUIRED) + add_library(dagmc-imported INTERFACE IMPORTED) link_directories(${DAGMC_LIBRARY_DIRS}) + target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES}) + target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS}) +else() + set(DAGMC_FOUND false) endif() #=============================================================================== # HDF5 for binary output #=============================================================================== -# Allow user to specify HDF5_ROOT -if (NOT (CMAKE_VERSION VERSION_LESS 3.12)) - cmake_policy(SET CMP0074 NEW) -endif() - # Unfortunately FindHDF5.cmake will always prefer a serial HDF5 installation # over a parallel installation if both appear on the user's PATH. To get around # this, we check for the environment variable HDF5_ROOT and if it exists, use it @@ -356,8 +361,7 @@ target_link_libraries(libopenmc ${ldflags} ${HDF5_LIBRARIES} ${HDF5_HL_LIBRARIES if(dagmc) target_compile_definitions(libopenmc PRIVATE DAGMC) - target_link_libraries(libopenmc ${DAGMC_LIBRARIES}) - target_include_directories(libopenmc PRIVATE ${DAGMC_INCLUDE_DIRS}) + target_link_libraries(libopenmc dagmc-imported) endif() #=============================================================================== @@ -387,6 +391,8 @@ add_custom_command(TARGET libopenmc POST_BUILD # Install executable, scripts, manpage, license #=============================================================================== +configure_file(cmake/OpenMCConfig.cmake.in "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" @ONLY) + set(INSTALL_CONFIGDIR ${CMAKE_INSTALL_LIBDIR}/cmake/OpenMC) install(TARGETS openmc libopenmc faddeeva EXPORT openmc-targets @@ -403,7 +409,7 @@ install(EXPORT openmc-targets # separately since they are managed by CMake) install(DIRECTORY vendor/faddeeva DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) install(DIRECTORY src/relaxng DESTINATION ${CMAKE_INSTALL_DATADIR}/openmc) -install(FILES cmake/OpenMCConfig.cmake DESTINATION ${INSTALL_CONFIGDIR}) +install(FILES "${CMAKE_BINARY_DIR}${CMAKE_FILES_DIRECTORY}/OpenMCConfig.cmake" DESTINATION ${INSTALL_CONFIGDIR}) install(FILES man/man1/openmc.1 DESTINATION ${CMAKE_INSTALL_MANDIR}/man1) install(FILES LICENSE DESTINATION "${CMAKE_INSTALL_DOCDIR}" RENAME copyright) install(DIRECTORY include/ DESTINATION ${CMAKE_INSTALL_INCLUDEDIR}) diff --git a/cmake/Modules/FindDAGMC.cmake b/cmake/Modules/FindDAGMC.cmake index c644c8886e..bb9bc56ddb 100644 --- a/cmake/Modules/FindDAGMC.cmake +++ b/cmake/Modules/FindDAGMC.cmake @@ -12,7 +12,9 @@ find_path(DAGMC_CMAKE_CONFIG NAMES DAGMCConfig.cmake PATHS ENV LD_LIBRARY_PATH PATH_SUFFIXES lib Lib cmake lib/cmake NO_DEFAULT_PATH) - -message(STATUS "Found DAGMC in ${DAGMC_CMAKE_CONFIG}") - -include(${DAGMC_CMAKE_CONFIG}/DAGMCConfig.cmake) +if(DAGMC_CMAKE_CONFIG) + message(STATUS "Found DAGMC in ${DAGMC_CMAKE_CONFIG}") + include(${DAGMC_CMAKE_CONFIG}/DAGMCConfig.cmake) +else() + message(WARNING "Cound not find DAGMC") +endif() diff --git a/cmake/OpenMCConfig.cmake b/cmake/OpenMCConfig.cmake.in similarity index 56% rename from cmake/OpenMCConfig.cmake rename to cmake/OpenMCConfig.cmake.in index 29a0e4542f..1c6a488f54 100644 --- a/cmake/OpenMCConfig.cmake +++ b/cmake/OpenMCConfig.cmake.in @@ -5,6 +5,15 @@ find_package(gsl-lite REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../gsl-lite) find_package(pugixml REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../pugixml) find_package(xtl REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtl) find_package(xtensor REQUIRED HINTS ${OpenMC_CMAKE_DIR}/../xtensor) +if(@DAGMC_FOUND@) + find_package(DAGMC REQUIRED HINTS @DAGMC_LIBRARY_DIRS@) + add_library(dagmc-imported INTERFACE IMPORTED) + target_link_libraries(dagmc-imported INTERFACE ${DAGMC_LIBRARIES}) + foreach(dir ${DAGMC_LIBRARY_DIRS}) + target_link_libraries(dagmc-imported INTERFACE "-L${dir}") + endforeach() + target_include_directories(dagmc-imported INTERFACE ${DAGMC_INCLUDE_DIRS}) +endif() if(NOT TARGET OpenMC::libopenmc) include("${OpenMC_CMAKE_DIR}/OpenMCTargets.cmake") diff --git a/openmc/data/photon.py b/openmc/data/photon.py index 32ed3010a4..9cc697b782 100644 --- a/openmc/data/photon.py +++ b/openmc/data/photon.py @@ -784,8 +784,9 @@ class IncidentPhoton(EqualityMixin): sub_group = shell_group.create_group(key) # Write atomic relaxation - if key in self.atomic_relaxation.subshells: - self.atomic_relaxation.to_hdf5(sub_group, key) + if self.atomic_relaxation is not None: + if key in self.atomic_relaxation.subshells: + self.atomic_relaxation.to_hdf5(sub_group, key) else: continue diff --git a/openmc/surface.py b/openmc/surface.py index 55e3f459d3..92c828d65f 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -1,7 +1,7 @@ from abc import ABCMeta, abstractmethod from collections import OrderedDict from copy import deepcopy -from numbers import Real, Integral +from numbers import Real from xml.etree import ElementTree as ET from warnings import warn @@ -19,6 +19,11 @@ _WARNING_UPPER = """\ will not accept the capitalized version.\ """ +_WARNING_KWARGS = """\ +"{}(...) accepts keyword arguments only for '{}'. Future versions of OpenMC \ +will not accept positional parameters for superclass arguments.\ +""" + class Surface(IDManagerMixin, metaclass=ABCMeta): """An implicit surface with an associated boundary condition. @@ -60,6 +65,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): next_id = 1 used_ids = set() + _atol = 1.e-12 def __init__(self, surface_id=None, boundary_type='transmission', name=''): self.id = surface_id @@ -179,13 +185,85 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): return memo[self] - @abstractmethod - def evaluate(self, point): - pass + def normalize(self, coeffs=None): + """Normalize coefficients by first nonzero value + + Parameters + ---------- + coeffs : tuple, optional + Tuple of surface coefficients to normalize. Defaults to None. If no + coefficients are supplied then the coefficients will be taken from + the current Surface. + + Returns + ------- + tuple of normalized coefficients + + """ + if coeffs is None: + coeffs = self._get_base_coeffs() + coeffs = np.asarray(coeffs) + nonzeros = ~np.isclose(coeffs, 0., rtol=0., atol=self._atol) + norm_factor = np.abs(coeffs[nonzeros][0]) + return tuple([c/norm_factor for c in coeffs]) + + def is_equal(self, other): + """Determine if this Surface is equivalent to another + + Parameters + ---------- + other : instance of openmc.Surface + Instance of openmc.Surface that should be compared to the current + surface + + """ + coeffs1 = self.normalize(self._get_base_coeffs()) + coeffs2 = self.normalize(other._get_base_coeffs()) + + return np.allclose(coeffs1, coeffs2, rtol=0., atol=self._atol) @abstractmethod - def translate(self, vector): - pass + def _get_base_coeffs(self): + """Return polynomial coefficients representing the implicit surface + equation. + + """ + + @abstractmethod + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + Evaluation of the surface polynomial at point :math:`(x',y',z')` + + """ + + @abstractmethod + def translate(self, vector, inplace=False): + """Translate surface in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which surface should be translated + inplace : boolean + Whether or not to return a new instance of this Surface or to + modify the coefficients of this Surface. Defaults to False + + Returns + ------- + instance of openmc.Surface + Translated surface + + """ def to_xml_element(self): """Return XML representation of the surface @@ -228,21 +306,7 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): # Determine appropriate class surf_type = elem.get('type') - surface_classes = { - 'plane': Plane, - 'x-plane': XPlane, - 'y-plane': YPlane, - 'z-plane': ZPlane, - 'x-cylinder': XCylinder, - 'y-cylinder': YCylinder, - 'z-cylinder': ZCylinder, - 'sphere': Sphere, - 'x-cone': XCone, - 'y-cone': YCone, - 'z-cone': ZCone, - 'quadric': Quadric, - } - cls = surface_classes[surf_type] + cls = _SURFACE_CLASSES[surf_type] # Determine ID, boundary type, coefficients kwargs = {} @@ -268,62 +332,102 @@ class Surface(IDManagerMixin, metaclass=ABCMeta): Instance of surface subclass """ + surface_id = int(group.name.split('/')[-1].lstrip('surface ')) name = group['name'][()].decode() if 'name' in group else '' surf_type = group['type'][()].decode() bc = group['boundary_type'][()].decode() coeffs = group['coefficients'][...] + kwargs = {'boundary_type': bc, 'name': name, 'surface_id': surface_id} - # Create the Surface based on its type - if surf_type == 'x-plane': - x0 = coeffs[0] - surface = XPlane(x0, bc, name, surface_id) + cls = _SURFACE_CLASSES[surf_type] - elif surf_type == 'y-plane': - y0 = coeffs[0] - surface = YPlane(y0, bc, name, surface_id) - - elif surf_type == 'z-plane': - z0 = coeffs[0] - surface = ZPlane(z0, bc, name, surface_id) - - elif surf_type == 'plane': - A, B, C, D = coeffs - surface = Plane(A, B, C, D, bc, name, surface_id) - - elif surf_type == 'x-cylinder': - y0, z0, r = coeffs - surface = XCylinder(y0, z0, r, bc, name, surface_id) - - elif surf_type == 'y-cylinder': - x0, z0, r = coeffs - surface = YCylinder(x0, z0, r, bc, name, surface_id) - - elif surf_type == 'z-cylinder': - x0, y0, r = coeffs - surface = ZCylinder(x0, y0, r, bc, name, surface_id) - - elif surf_type == 'sphere': - x0, y0, z0, r = coeffs - surface = Sphere(x0, y0, z0, r, bc, name, surface_id) - - elif surf_type in ['x-cone', 'y-cone', 'z-cone']: - x0, y0, z0, r2 = coeffs - if surf_type == 'x-cone': - surface = XCone(x0, y0, z0, r2, bc, name, surface_id) - elif surf_type == 'y-cone': - surface = YCone(x0, y0, z0, r2, bc, name, surface_id) - elif surf_type == 'z-cone': - surface = ZCone(x0, y0, z0, r2, bc, name, surface_id) - - elif surf_type == 'quadric': - a, b, c, d, e, f, g, h, j, k = coeffs - surface = Quadric(a, b, c, d, e, f, g, h, j, k, bc, name, surface_id) - - return surface + return cls(*coeffs, **kwargs) -class Plane(Surface): +class PlaneMixin(metaclass=ABCMeta): + """A Plane mixin class for all operations on order 1 surfaces""" + def __init__(self, **kwargs): + super().__init__(**kwargs) + self._periodic_surface = None + + @property + def periodic_surface(self): + return self._periodic_surface + + @periodic_surface.setter + def periodic_surface(self, periodic_surface): + check_type('periodic surface', periodic_surface, Plane) + self._periodic_surface = periodic_surface + periodic_surface._periodic_surface = self + + def evaluate(self, point): + """Evaluate the surface equation at a given point. + + Parameters + ---------- + point : 3-tuple of float + The Cartesian coordinates, :math:`(x',y',z')`, at which the surface + equation should be evaluated. + + Returns + ------- + float + :math:`Ax' + By' + Cz' - D` + + """ + + x, y, z = point + a, b, c, d = self._get_base_coeffs() + return a*x + b*y + c*z - d + + def translate(self, vector, inplace=False): + """Translate surface in given direction + + Parameters + ---------- + vector : iterable of float + Direction in which surface should be translated + inplace : boolean + Whether or not to return a new instance of a Plane or to modify the + coefficients of this plane. Defaults to False + + Returns + ------- + openmc.Plane + Translated surface + + """ + vx, vy, vz = vector + a, b, c, d = self._get_base_coeffs() + d = d + a*vx + b*vy + c*vz + + surf = self if inplace else self.clone() + + setattr(surf, surf._coeff_keys[-1], d) + + return surf + + def to_xml_element(self): + """Return XML representation of the surface + + Returns + ------- + element : xml.etree.ElementTree.Element + XML element containing source data + + """ + element = super().to_xml_element() + + # Add periodic surface pair information + if self.boundary_type == 'periodic': + if self.periodic_surface is not None: + element.set("periodic_surface_id", + str(self.periodic_surface.id)) + return element + + +class Plane(PlaneMixin, Surface): """An arbitrary plane of the form :math:`Ax + By + Cz = D`. Parameters @@ -376,19 +480,34 @@ class Plane(Surface): _type = 'plane' _coeff_keys = ('a', 'b', 'c', 'd') - def __init__(self, a=1., b=0., c=0., d=0., boundary_type='transmission', - name='', surface_id=None, **kwargs): - super().__init__(surface_id, boundary_type, name=name) - self._periodic_surface = None - self.a = a - self.b = b - self.c = c - self.d = d - for k, v in kwargs.items(): - if k in 'ABCD': - warn(_WARNING_UPPER.format(type(self).__name__, k.lower(), k), - FutureWarning) - setattr(self, k.lower(), v) + def __init__(self, a=1., b=0., c=0., d=0., *args, **kwargs): + # *args should ultimately be limited to a, b, c, d as specified in + # __init__, but to preserve the API it is allowed to accept Surface + # parameters for now, but will raise warnings if this is done. + argtup = ('a', 'b', 'c', 'd', 'boundary_type', 'name', 'surface_id') + kwargs.update(dict(zip(argtup, args))) + + # Warn if Surface parameters are passed by position, not by keyword + superkwargs = {} + for k in ('boundary_type', 'name', 'surface_id'): + val = kwargs.get(k, None) + if val is not None: + superkwargs[k] = val + warn(_WARNING_KWARGS.format(type(self), k), + FutureWarning) + + super().__init__(**superkwargs) + + for key, val in zip(self._coeff_keys, (a, b, c, d)): + setattr(self, key, val) + + # Warn if capital letter arguments are passed + for k in 'ABCD': + val = kwargs.pop(k, None) + if val is not None: + warn(_WARNING_UPPER.format(type(self), k.lower(), k), + FutureWarning) + setattr(self, k.lower(), val) @property def a(self): @@ -406,10 +525,6 @@ class Plane(Surface): def d(self): return self.coefficients['d'] - @property - def periodic_surface(self): - return self._periodic_surface - @a.setter def a(self, a): check_type('A coefficient', a, Real) @@ -430,68 +545,8 @@ class Plane(Surface): check_type('D coefficient', d, Real) self._coefficients['d'] = d - @periodic_surface.setter - def periodic_surface(self, periodic_surface): - check_type('periodic surface', periodic_surface, Plane) - self._periodic_surface = periodic_surface - periodic_surface._periodic_surface = self - - def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`Ax' + By' + Cz' - D` - - """ - - x, y, z = point - return self.a*x + self.b*y + self.c*z - self.d - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.Plane - Translated surface - - """ - vx, vy, vz = vector - d = self.d + self.a*vx + self.b*vy + self.c*vz - if d == self.d: - return self - else: - return type(self)(a=self.a, b=self.b, c=self.c, d=d) - - def to_xml_element(self): - """Return XML representation of the surface - - Returns - ------- - element : xml.etree.ElementTree.Element - XML element containing source data - - """ - element = super().to_xml_element() - - # Add periodic surface pair information - if self.boundary_type == 'periodic': - if self.periodic_surface is not None: - element.set("periodic_surface_id", str(self.periodic_surface.id)) - return element + def _get_base_coeffs(self): + return (self.a, self.b, self.c, self.d) @classmethod def from_points(cls, p1, p2, p3, **kwargs): @@ -526,7 +581,7 @@ class Plane(Surface): return cls(a=a, b=b, c=c, d=d, **kwargs) -class XPlane(Plane): +class XPlane(PlaneMixin, Surface): """A plane perpendicular to the x axis of the form :math:`x - x_0 = 0` Parameters @@ -568,9 +623,15 @@ class XPlane(Plane): _type = 'x-plane' _coeff_keys = ('x0',) - def __init__(self, x0=0., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id=surface_id, boundary_type=boundary_type, name=name) + def __init__(self, x0=0., *args, **kwargs): + # work around for accepting Surface kwargs as positional parameters + # until they are deprecated + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) + for k in argsdict: + warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + kwargs.update(argsdict) + + super().__init__(**kwargs) self.x0 = x0 @property @@ -582,30 +643,10 @@ class XPlane(Plane): check_type('x0 coefficient', x0, Real) self._coefficients['x0'] = x0 + def _get_base_coeffs(self): + return (1., 0., 0., self.x0) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the x-plane surface, the - half-spaces are unbounded in their y- and z- directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([-np.inf, -np.inf, -np.inf]), np.array([self.x0, np.inf, np.inf])) @@ -614,44 +655,13 @@ class XPlane(Plane): np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`x' - x_0` - - """ return point[0] - self.x0 - def translate(self, vector): - """Translate surface in given direction - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.XPlane - Translated surface - - """ - vx = vector[0] - if vx == 0: - return self - else: - return type(self)(x0=self.x0 + vx) +Plane.register(XPlane) -class YPlane(Plane): +class YPlane(PlaneMixin, Surface): """A plane perpendicular to the y axis of the form :math:`y - y_0 = 0` Parameters @@ -693,9 +703,15 @@ class YPlane(Plane): _type = 'y-plane' _coeff_keys = ('y0',) - def __init__(self, y0=0., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id=surface_id, boundary_type=boundary_type, name=name) + def __init__(self, y0=0., *args, **kwargs): + # work around for accepting Surface kwargs as positional parameters + # until they are deprecated + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) + for k in argsdict: + warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + kwargs.update(argsdict) + + super().__init__(**kwargs) self.y0 = y0 @property @@ -707,30 +723,10 @@ class YPlane(Plane): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 + def _get_base_coeffs(self): + return (0., 1., 0., self.y0) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the y-plane surface, the - half-spaces are unbounded in their x- and z- directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, self.y0, np.inf])) @@ -739,60 +735,29 @@ class YPlane(Plane): np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`y' - y_0` - - """ return point[1] - self.y0 - def translate(self, vector): - """Translate surface in given direction - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.YPlane - Translated surface - - """ - vy = vector[1] - if vy == 0.0: - return self - else: - return type(self)(y0=self.y0 + vy) +Plane.register(YPlane) -class ZPlane(Plane): +class ZPlane(PlaneMixin, Surface): """A plane perpendicular to the z axis of the form :math:`z - z_0 = 0` Parameters ---------- - surface_id : int, optional - Unique identifier for the surface. If not specified, an identifier will - automatically be assigned. + z0 : float, optional + Location of the plane. Defaults to 0. boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. Only axis-aligned periodicity is supported, i.e., x-planes can only be paired with x-planes. - z0 : float, optional - Location of the plane. Defaults to 0. name : str, optional Name of the plane. If not specified, the name will be the empty string. + surface_id : int, optional + Unique identifier for the surface. If not specified, an identifier will + automatically be assigned. Attributes ---------- @@ -818,9 +783,15 @@ class ZPlane(Plane): _type = 'z-plane' _coeff_keys = ('z0',) - def __init__(self, z0=0., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id=surface_id, boundary_type=boundary_type, name=name) + def __init__(self, z0=0., *args, **kwargs): + # work around for accepting Surface kwargs as positional parameters + # until they are deprecated + argsdict = dict(zip(('boundary_type', 'name', 'surface_id'), args)) + for k in argsdict: + warn(_WARNING_KWARGS.format(type(self).__name__, k), FutureWarning) + kwargs.update(argsdict) + + super().__init__(**kwargs) self.z0 = z0 @property @@ -832,30 +803,10 @@ class ZPlane(Plane): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + def _get_base_coeffs(self): + return (0., 0., 1., self.z0) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the z-plane surface, the - half-spaces are unbounded in their x- and y- directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([-np.inf, -np.inf, -np.inf]), np.array([np.inf, np.inf, self.z0])) @@ -863,6 +814,58 @@ class ZPlane(Plane): return (np.array([-np.inf, -np.inf, self.z0]), np.array([np.inf, np.inf, np.inf])) + def evaluate(self, point): + return point[2] - self.z0 + + +Plane.register(ZPlane) + + +class QuadricMixin(metaclass=ABCMeta): + """A Mixin class implementing common functionality for quadric surfaces""" + + def get_Abc(self, coeffs=None): + """Compute matrix, vector, and scalar coefficients for this surface or + for a specified set of coefficients. + + Parameters + ---------- + coeffs : tuple, optional + Tuple of coefficients from which to compute the quadric elements. + If none are supplied the coefficients of this surface will be used. + """ + if coeffs is None: + a, b, c, d, e, f, g, h, j, k = self._get_base_coeffs() + else: + a, b, c, d, e, f, g, h, j, k = coeffs + + A = np.array([[a, d/2, f/2], [d/2, b, e/2], [f/2, e/2, c]]) + bvec = np.array([g, h, j]) + + return A, bvec, k + + def eigh(self, coeffs=None): + """Wrapper method for returning eigenvalues and eigenvectors of this + quadric surface which is used for transformations. + + Parameters + ---------- + coeffs : tuple, optional + Tuple of coefficients from which to compute the quadric elements. + If none are supplied the coefficients of this surface will be used. + + Returns + ------- + w, v : tuple of numpy arrays with shapes (3,) and (3,3) respectively + Returns the eigenvalues and eigenvectors of the quadric matrix A + that represents the supplied coefficients. The vector w contains + the eigenvalues in ascending order and the matrix v contains the + eigenvectors such that v[:,i] is the eigenvector corresponding to + the eigenvalue w[i]. + + """ + return np.linalg.eigh(self.get_Abc(coeffs=coeffs)[0]) + def evaluate(self, point): """Evaluate the surface equation at a given point. @@ -875,37 +878,269 @@ class ZPlane(Plane): Returns ------- float - :math:`z' - z_0` + :math:`Ax'^2 + By'^2 + Cz'^2 + Dx'y' + Ey'z' + Fx'z' + Gx' + Hy' + + Jz' + K = 0` """ - return point[2] - self.z0 + x = np.asarray(point) + A, b, c = self.get_Abc() + return np.matmul(x.T, np.matmul(A, x)) + np.matmul(b.T, x) + c - def translate(self, vector): + def translate(self, vector, inplace=False): """Translate surface in given direction Parameters ---------- vector : iterable of float Direction in which surface should be translated + inplace : boolean + Whether to return a clone of the Surface or the Surface itself. + Defaults to False Returns ------- - openmc.ZPlane + openmc.Surface Translated surface """ - vz = vector[2] - if vz == 0.0: - return self + vector = np.asarray(vector) + + surf = self if inplace else self.clone() + + if set(('x0', 'y0', 'z0')).intersection(set(surf._coeff_keys)): + for vi, xi in zip(vector, ('x0', 'y0', 'z0')): + val = getattr(surf, xi, None) + if val is not None: + setattr(surf, xi, val + vi) else: - return type(self)(z0=self.z0 + vz) + A, bvec, cnst = self.get_Abc() + + g, h, j = bvec - 2*np.matmul(vector.T, A) + k = cnst + np.matmul(vector.T, np.matmul(A, vector)) \ + - np.matmul(bvec.T, vector) + + for key, val in zip(('g', 'h', 'j', 'k'), (g, h, j, k)): + setattr(surf, key, val) + + return surf -class Cylinder(Surface): - """A cylinder whose length is parallel to the x-, y-, or z-axis. +class Cylinder(QuadricMixin, Surface): + """A cylinder with radius r, centered on the point (x0, y0, z0) with an + axis specified by the line through points (x0, y0, z0) and (x0+dx, y0+dy, + z0+dz) Parameters ---------- + x0 : float, optional + x-coordinate for the origin of the Cylinder. Defaults to 0 + y0 : float, optional + y-coordinate for the origin of the Cylinder. Defaults to 0 + z0 : float, optional + z-coordinate for the origin of the Cylinder. Defaults to 0 + r : float, optional + Radius of the cylinder. Defaults to 1. + dx : float, optional + x-component of the vector representing the axis of the cylinder. + Defaults to 0. + dy : float, optional + y-component of the vector representing the axis of the cylinder. + Defaults to 0. + dz : float, optional + z-component of the vector representing the axis of the cylinder. + Defaults to 1. + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional + Boundary condition that defines the behavior for particles hitting the + surface. Defaults to transmissive boundary condition where particles + freely pass through the surface. + name : str, optional + Name of the cylinder. If not specified, the name will be the empty + string. + surface_id : int, optional + Unique identifier for the surface. If not specified, an identifier will + automatically be assigned. + + Attributes + ---------- + x0 : float + x-coordinate for the origin of the Cylinder + y0 : float + y-coordinate for the origin of the Cylinder + z0 : float + z-coordinate for the origin of the Cylinder + r : float + Radius of the cylinder + dx : float + x-component of the vector representing the axis of the cylinder + dy : float + y-component of the vector representing the axis of the cylinder + dz : float + z-component of the vector representing the axis of the cylinder + boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} + Boundary condition that defines the behavior for particles hitting the + surface. + coefficients : dict + Dictionary of surface coefficients + id : int + Unique identifier for the surface + name : str + Name of the surface + type : str + Type of the surface + + """ + _type = 'cylinder' + _coeff_keys = ('x0', 'y0', 'z0', 'r', 'dx', 'dy', 'dz') + + def __init__(self, x0=0., y0=0., z0=0., r=1., dx=0., dy=0., dz=1., **kwargs): + raise NotImplementedError('There is no C++ implementation for general ' + 'Cylinders yet, please use ' + 'openmc.model.funcs.cylinder_from_points to ' + 'return a Quadric instance instead for now') + + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r, dx, dy, dz)): + setattr(self, key, val) + + @property + def x0(self): + return self.coefficients['x0'] + + @property + def y0(self): + return self.coefficients['y0'] + + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r(self): + return self.coefficients['r'] + + @property + def dx(self): + return self.coefficients['dx'] + + @property + def dy(self): + return self.coefficients['dy'] + + @property + def dz(self): + return self.coefficients['dz'] + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + + @dx.setter + def dx(self, dx): + check_type('dx coefficient', dx, Real) + self._coefficients['dx'] = dx + + @dy.setter + def dy(self, dy): + check_type('dy coefficient', dy, Real) + self._coefficients['dy'] = dy + + @dz.setter + def dz(self, dz): + check_type('dz coefficient', dz, Real) + self._coefficients['dz'] = dz + + def _get_base_coeffs(self): + # Get x, y, z coordinates of two points + x1, y1, z1 = self.x0, self.y0, self.z0 + x2, y2, z2 = x1 + self.dx, y1 + self.dy, z1 + self.dz + r = self.r + + # Define intermediate terms + dx = x2 - x1 + dy = y2 - y1 + dz = z2 - z1 + cx = y1*z2 - y2*z1 + cy = x2*z1 - x1*z2 + cz = x1*y2 - x2*y1 + + # Given p=(x,y,z), p1=(x1, y1, z1), p2=(x2, y2, z2), the equation + # for the cylinder can be derived as + # r = |(p - p1) тип (p - p2)| / |p2 - p1|. + # Expanding out all terms and grouping according to what Quadric + # expects gives the following coefficients. + a = dy*dy + dz*dz + b = dx*dx + dz*dz + c = dx*dx + dy*dy + d = -2*dx*dy + e = -2*dy*dz + f = -2*dx*dz + g = 2*(cy*dz - cz*dy) + h = 2*(cz*dx - cx*dz) + j = 2*(cx*dy - cy*dx) + k = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r + + return (a, b, c, d, e, f, g, h, j, k) + + @classmethod + def from_points(cls, p1, p2, r=1., **kwargs): + """Return a cylinder given points that define the axis and a radius. + + Parameters + ---------- + p1, p2 : 3-tuples + Points that pass through the plane, p1 will be used as (x0, y0, z0) + r : float, optional + Radius of the cylinder. Defaults to 1. + kwargs : dict + Keyword arguments passed to the :class:`Cylinder` constructor + + Returns + ------- + Cylinder + Cylinder that has an axis through the points p1 and p2, and a + radius r. + + """ + raise NotImplementedError('There is no C++ implementation for general ' + 'Cylinders yet, please use ' + 'openmc.model.funcs.cylinder_from_points to ' + 'return a Quadric instance instead for now') + # Convert to numpy arrays + p1 = np.asarray(p1) + p2 = np.asarray(p2) + x0, y0, z0 = p1 + dx, dy, dz = p2 - p1 + + return cls(x0=x0, y0=y0, z0=z0, r=r, dx=dx, dy=dy, dz=dz, **kwargs) + + +class XCylinder(QuadricMixin, Surface): + """An infinite cylinder whose length is parallel to the x-axis of the form + :math:`(y - y_0)^2 + (z - z_0)^2 = r^2`. + + Parameters + ---------- + y0 : float, optional + y-coordinate for the origin of the Cylinder. Defaults to 0 + z0 : float, optional + z-coordinate for the origin of the Cylinder. Defaults to 0 r : float, optional Radius of the cylinder. Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional @@ -921,6 +1156,10 @@ class Cylinder(Surface): Attributes ---------- + y0 : float + y-coordinate for the origin of the Cylinder + z0 : float + z-coordinate for the origin of the Cylinder r : float Radius of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} @@ -935,76 +1174,21 @@ class Cylinder(Surface): type : str Type of the surface - """ - def __init__(self, r=1., boundary_type='transmission', - name='', surface_id=None): - super().__init__(surface_id, boundary_type, name=name) - self.r = r - - @property - def r(self): - return self.coefficients['r'] - - @r.setter - def r(self, r): - check_type('r coefficient', r, Real) - self._coefficients['r'] = r - - -class XCylinder(Cylinder): - """An infinite cylinder whose length is parallel to the x-axis of the form - :math:`(y - y_0)^2 + (z - z_0)^2 = r^2`. - - Parameters - ---------- - y0 : float, optional - y-coordinate of the center of the cylinder. Defaults to 0. - z0 : float, optional - z-coordinate of the center of the cylinder. Defaults to 0. - r : float, optional - Radius of the cylinder. Defaults to 0. - boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional - Boundary condition that defines the behavior for particles hitting the - surface. Defaults to transmissive boundary condition where particles - freely pass through the surface. - name : str, optional - Name of the cylinder. If not specified, the name will be the empty - string. - surface_id : int, optional - Unique identifier for the surface. If not specified, an identifier will - automatically be assigned. - - Attributes - ---------- - y0 : float - y-coordinate of the center of the cylinder - z0 : float - z-coordinate of the center of the cylinder - boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} - Boundary condition that defines the behavior for particles hitting the - surface. - coefficients : dict - Dictionary of surface coefficients - id : int - Unique identifier for the surface - name : str - Name of the surface - type : str - Type of the surface - """ _type = 'x-cylinder' _coeff_keys = ('y0', 'z0', 'r') - def __init__(self, y0=0., z0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, y0=0., z0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R - super().__init__(r, boundary_type, name, surface_id) - self.y0 = y0 - self.z0 = z0 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (y0, z0, r)): + setattr(self, key, val) @property def y0(self): @@ -1014,6 +1198,10 @@ class XCylinder(Cylinder): def z0(self): return self.coefficients['z0'] + @property + def r(self): + return self.coefficients['r'] + @y0.setter def y0(self, y0): check_type('y0 coefficient', y0, Real) @@ -1024,31 +1212,21 @@ class XCylinder(Cylinder): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + + def _get_base_coeffs(self): + y0, z0, r = self.y0, self.z0, self.r + + a = d = e = f = g = 0. + b = c = 1. + h, j, k = -2*y0, -2*z0, y0*y0 + z0*z0 - r*r + + return (a, b, c, d, e, f, g, h, j, k) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the x-cylinder surface, - the negative half-space is unbounded in the x- direction and the - positive half-space is unbounded in all directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([-np.inf, self.y0 - self.r, self.z0 - self.r]), np.array([np.inf, self.y0 + self.r, self.z0 + self.r])) @@ -1057,57 +1235,24 @@ class XCylinder(Cylinder): np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(y' - y_0)^2 + (z' - z_0)^2 - r^2` - - """ y = point[1] - self.y0 z = point[2] - self.z0 - return y**2 + z**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.XCylinder - Translated surface - - """ - vx, vy, vz = vector - if vy == 0.0 and vz == 0.0: - return self - else: - y0 = self.y0 + vy - z0 = self.z0 + vz - return type(self)(y0=y0, z0=z0, r=self.r) + return y*y + z*z - self.r**2 -class YCylinder(Cylinder): +Cylinder.register(XCylinder) + + +class YCylinder(QuadricMixin, Surface): """An infinite cylinder whose length is parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 = r^2`. Parameters ---------- x0 : float, optional - x-coordinate of the center of the cylinder. Defaults to 0. + x-coordinate for the origin of the Cylinder. Defaults to 0 z0 : float, optional - z-coordinate of the center of the cylinder. Defaults to 0. + z-coordinate for the origin of the Cylinder. Defaults to 0 r : float, optional Radius of the cylinder. Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional @@ -1124,9 +1269,11 @@ class YCylinder(Cylinder): Attributes ---------- x0 : float - x-coordinate of the center of the cylinder + x-coordinate for the origin of the Cylinder z0 : float - z-coordinate of the center of the cylinder + z-coordinate for the origin of the Cylinder + r : float + Radius of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -1144,14 +1291,16 @@ class YCylinder(Cylinder): _type = 'y-cylinder' _coeff_keys = ('x0', 'z0', 'r') - def __init__(self, x0=0., z0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, x0=0., z0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R - super().__init__(r, boundary_type, name, surface_id) - self.x0 = x0 - self.z0 = z0 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, z0, r)): + setattr(self, key, val) @property def x0(self): @@ -1161,6 +1310,10 @@ class YCylinder(Cylinder): def z0(self): return self.coefficients['z0'] + @property + def r(self): + return self.coefficients['r'] + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1171,31 +1324,21 @@ class YCylinder(Cylinder): check_type('z0 coefficient', z0, Real) self._coefficients['z0'] = z0 + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + + def _get_base_coeffs(self): + x0, z0, r = self.x0, self.z0, self.r + + b = d = e = f = h = 0. + a = c = 1. + g, j, k = -2*x0, -2*z0, x0*x0 + z0*z0 - r*r + + return (a, b, c, d, e, f, g, h, j, k) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the y-cylinder surface, - the negative half-space is unbounded in the y- direction and the - positive half-space is unbounded in all directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([self.x0 - self.r, -np.inf, self.z0 - self.r]), np.array([self.x0 + self.r, np.inf, self.z0 + self.r])) @@ -1204,76 +1347,45 @@ class YCylinder(Cylinder): np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (z' - z_0)^2 - r^2` - - """ x = point[0] - self.x0 z = point[2] - self.z0 - return x**2 + z**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.YCylinder - Translated surface - - """ - vx, vy, vz = vector - if vx == 0.0 and vz == 0.0: - return self - else: - x0 = self.x0 + vx - z0 = self.z0 + vz - return type(self)(x0=x0, z0=z0, r=self.r) + return x*x + z*z - self.r**2 -class ZCylinder(Cylinder): +Cylinder.register(YCylinder) + + +class ZCylinder(QuadricMixin, Surface): """An infinite cylinder whose length is parallel to the z-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 = r^2`. Parameters ---------- - surface_id : int, optional - Unique identifier for the surface. If not specified, an identifier will - automatically be assigned. + x0 : float, optional + x-coordinate for the origin of the Cylinder. Defaults to 0 + y0 : float, optional + y-coordinate for the origin of the Cylinder. Defaults to 0 + r : float, optional + Radius of the cylinder. Defaults to 1. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'}, optional Boundary condition that defines the behavior for particles hitting the surface. Defaults to transmissive boundary condition where particles freely pass through the surface. - x0 : float, optional - x-coordinate of the center of the cylinder. Defaults to 0. - y0 : float, optional - y-coordinate of the center of the cylinder. Defaults to 0. - r : float, optional - Radius of the cylinder. Defaults to 1. name : str, optional Name of the cylinder. If not specified, the name will be the empty string. + surface_id : int, optional + Unique identifier for the surface. If not specified, an identifier will + automatically be assigned. Attributes ---------- x0 : float - x-coordinate of the center of the cylinder + x-coordinate for the origin of the Cylinder y0 : float - y-coordinate of the center of the cylinder + y-coordinate for the origin of the Cylinder + r : float + Radius of the cylinder boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -1291,14 +1403,16 @@ class ZCylinder(Cylinder): _type = 'z-cylinder' _coeff_keys = ('x0', 'y0', 'r') - def __init__(self, x0=0., y0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, x0=0., y0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R - super().__init__(r, boundary_type, name, surface_id) - self.x0 = x0 - self.y0 = y0 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, r)): + setattr(self, key, val) @property def x0(self): @@ -1308,6 +1422,10 @@ class ZCylinder(Cylinder): def y0(self): return self.coefficients['y0'] + @property + def r(self): + return self.coefficients['r'] + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1318,31 +1436,21 @@ class ZCylinder(Cylinder): check_type('y0 coefficient', y0, Real) self._coefficients['y0'] = y0 + @r.setter + def r(self, r): + check_type('r coefficient', r, Real) + self._coefficients['r'] = r + + def _get_base_coeffs(self): + x0, y0, r = self.x0, self.y0, self.r + + c = d = e = f = j = 0. + a = b = 1. + g, h, k = -2*x0, -2*y0, x0*x0 + y0*y0 - r*r + + return (a, b, c, d, e, f, g, h, j, k) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. For the z-cylinder surface, - the negative half-space is unbounded in the z- direction and the - positive half-space is unbounded in all directions. To represent - infinity, numpy.inf is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([self.x0 - self.r, self.y0 - self.r, -np.inf]), np.array([self.x0 + self.r, self.y0 + self.r, np.inf])) @@ -1351,48 +1459,15 @@ class ZCylinder(Cylinder): np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (y' - y_0)^2 - r^2` - - """ x = point[0] - self.x0 y = point[1] - self.y0 - return x**2 + y**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.ZCylinder - Translated surface - - """ - vx, vy, vz = vector - if vx == 0.0 and vy == 0.0: - return self - else: - x0 = self.x0 + vx - y0 = self.y0 + vy - return type(self)(x0=x0, y0=y0, r=self.r) + return x*x + y*y - self.r**2 -class Sphere(Surface): +Cylinder.register(ZCylinder) + + +class Sphere(QuadricMixin, Surface): """A sphere of the form :math:`(x - x_0)^2 + (y - y_0)^2 + (z - z_0)^2 = r^2`. Parameters @@ -1442,16 +1517,16 @@ class Sphere(Surface): _type = 'sphere' _coeff_keys = ('x0', 'y0', 'z0', 'r') - def __init__(self, x0=0., y0=0., z0=0., r=1., boundary_type='transmission', - name='', surface_id=None, *, R=None): + def __init__(self, x0=0., y0=0., z0=0., r=1., **kwargs): + R = kwargs.pop('R', None) if R is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r', 'R'), + FutureWarning) r = R - super().__init__(surface_id, boundary_type, name=name) - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.r = r + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r)): + setattr(self, key, val) @property def x0(self): @@ -1489,30 +1564,16 @@ class Sphere(Surface): check_type('r coefficient', r, Real) self._coefficients['r'] = r + def _get_base_coeffs(self): + x0, y0, z0, r = self.x0, self.y0, self.z0, self.r + a = b = c = 1. + d = e = f = 0. + g, h, j = -2*x0, -2*y0, -2*z0 + k = x0*x0 + y0*y0 + z0*z0 - r*r + + return (a, b, c, d, e, f, g, h, j, k) + def bounding_box(self, side): - """Determine an axis-aligned bounding box. - - An axis-aligned bounding box for surface half-spaces is represented by - its lower-left and upper-right coordinates. The positive half-space of a - sphere is unbounded in all directions. To represent infinity, numpy.inf - is used. - - Parameters - ---------- - side : {'+', '-'} - Indicates the negative or positive half-space - - Returns - ------- - numpy.ndarray - Lower-left coordinates of the axis-aligned bounding box for the - desired half-space - numpy.ndarray - Upper-right coordinates of the axis-aligned bounding box for the - desired half-space - - """ - if side == '-': return (np.array([self.x0 - self.r, self.y0 - self.r, self.z0 - self.r]), @@ -1523,50 +1584,13 @@ class Sphere(Surface): np.array([np.inf, np.inf, np.inf])) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (y' - y_0)^2 + (z' - z_0)^2 - r^2` - - """ x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return x**2 + y**2 + z**2 - self.r**2 - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.Sphere - Translated surface - - """ - vx, vy, vz = vector - if vx == 0.0 and vy == 0.0 and vz == 0.0: - return self - else: - x0 = self.x0 + vx - y0 = self.y0 + vy - z0 = self.z0 + vz - return type(self)(x0=x0, y0=y0, z0=z0, r=self.r) + return x*x + y*y + z*z - self.r**2 -class Cone(Surface): +class Cone(QuadricMixin, Surface): """A conical surface parallel to the x-, y-, or z-axis. Parameters @@ -1579,6 +1603,15 @@ class Cone(Surface): z-coordinate of the apex. Defaults to 0. r2 : float, optional Parameter related to the aperature. Defaults to 1. + dx : float, optional + x-component of the vector representing the axis of the cone. + Defaults to 0. + dy : float, optional + y-component of the vector representing the axis of the cone. + Defaults to 0. + dz : float, optional + z-component of the vector representing the axis of the cone. + Defaults to 1. surface_id : int, optional Unique identifier for the surface. If not specified, an identifier will automatically be assigned. @@ -1599,6 +1632,12 @@ class Cone(Surface): z-coordinate of the apex r2 : float Parameter related to the aperature + dx : float + x-component of the vector representing the axis of the cone. + dy : float + y-component of the vector representing the axis of the cone. + dz : float + z-component of the vector representing the axis of the cone. boundary_type : {'transmission, 'vacuum', 'reflective', 'white'} Boundary condition that defines the behavior for particles hitting the surface. @@ -1613,18 +1652,22 @@ class Cone(Surface): """ - _coeff_keys = ('x0', 'y0', 'z0', 'r2') + _type = 'cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2', 'dx', 'dy', 'dz') - def __init__(self, x0=0., y0=0., z0=0., r2=1., boundary_type='transmission', - name='', surface_id=None, *, R2=None): + def __init__(self, x0=0., y0=0., z0=0., r2=1., dx=0., dy=0., dz=1., **kwargs): + raise NotImplementedError('There is no C++ implementation for general ' + 'Cones yet, this functionality should be ' + 'added soon.') + R2 = kwargs.pop('R2', None) if R2 is not None: - warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), FutureWarning) + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) r2 = R2 - super().__init__(surface_id, boundary_type, name=name) - self.x0 = x0 - self.y0 = y0 - self.z0 = z0 - self.r2 = r2 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2, dx, dy, dz)): + setattr(self, key, val) @property def x0(self): @@ -1642,6 +1685,18 @@ class Cone(Surface): def r2(self): return self.coefficients['r2'] + @property + def dx(self): + return self.coefficients['dx'] + + @property + def dy(self): + return self.coefficients['dy'] + + @property + def dz(self): + return self.coefficients['dz'] + @x0.setter def x0(self, x0): check_type('x0 coefficient', x0, Real) @@ -1662,31 +1717,57 @@ class Cone(Surface): check_type('r^2 coefficient', r2, Real) self._coefficients['r2'] = r2 - def translate(self, vector): - """Translate surface in given direction + @dx.setter + def dx(self, dx): + check_type('dx coefficient', dx, Real) + self._coefficients['dx'] = dx - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated + @dy.setter + def dy(self, dy): + check_type('dy coefficient', dy, Real) + self._coefficients['dy'] = dy - Returns - ------- - openmc.Cone - Translated surface + @dz.setter + def dz(self, dz): + check_type('dz coefficient', dz, Real) + self._coefficients['dz'] = dz - """ - vx, vy, vz = vector - if vx == 0.0 and vy == 0.0 and vz == 0.0: - return self - else: - x0 = self.x0 + vx - y0 = self.y0 + vy - z0 = self.z0 + vz - return type(self)(x0=x0, y0=y0, z0=z0, r2=self.r2) + def _get_base_coeffs(self): + # The equation for a general cone with vertex at point p = (x0, y0, z0) + # and axis specified by the unit vector d = (dx, dy, dz) and opening + # half angle theta can be described by the equation + # + # (d*(r - p))^2 - (r - p)*(r - p)cos^2(theta) = 0 + # + # where * is the dot product and the vector r is the evaulation point + # r = (x, y, z) + # + # The argument r2 for cones is actually tan^2(theta) so that + # cos^2(theta) = 1 / (1 + r2) + + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + dx, dy, dz = self.dx, self.dy, self.dz + dnorm = dx*dx + dy*dy + dz*dz + dx /= dnorm + dy /= dnorm + dz /= dnorm + cos2 = 1 / (1 + r2) + + a = dx*dx - cos2 + b = dy*dy - cos2 + c = dz*dz - cos2 + d = 2*dx*dy + e = 2*dy*dz + f = 2*dx*dz + g = -2*(dx*dx*x0 + dx*dy*y0 + dx*dz*z0 - cos2) + h = -2*(dy*dy*y0 + dx*dy*x0 + dy*dz*z0 - cos2) + j = -2*(dz*dz*y0 + dx*dz*x0 + dy*dz*y0 - cos2) + k = (dx*x0 + dy*y0 + dz*z0)**2 - cos2*(x0*x0 + y0*y0 + z0*z0) + + return (a, b, c, d, e, f, g, h, j, k) -class XCone(Cone): +class XCone(QuadricMixin, Surface): """A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 = r^2 (x - x_0)^2`. @@ -1735,29 +1816,77 @@ class XCone(Cone): """ _type = 'x-cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2') + + def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + R2 = kwargs.pop('R2', None) + if R2 is not None: + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) + r2 = R2 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) + + @property + def x0(self): + return self.coefficients['x0'] + + @property + def y0(self): + return self.coefficients['y0'] + + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r2(self): + return self.coefficients['r2'] + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r2.setter + def r2(self, r2): + check_type('r^2 coefficient', r2, Real) + self._coefficients['r2'] = r2 + + def _get_base_coeffs(self): + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + + a = -r2 + b = c = 1. + d = e = f = 0. + g, h, j = 2*x0*r2, -2*y0, -2*z0 + k = y0*y0 + z0*z0 - r2*x0*x0 + + return (a, b, c, d, e, f, g, h, j, k) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(y' - y_0)^2 + (z' - z_0)^2 - r^2(x' - x_0)^2` - - """ x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return y**2 + z**2 - self.r2*x**2 + return y*y + z*z - self.r2*x*x -class YCone(Cone): +Cone.register(XCone) + + +class YCone(QuadricMixin, Surface): """A cone parallel to the y-axis of the form :math:`(x - x_0)^2 + (z - z_0)^2 = r^2 (y - y_0)^2`. @@ -1806,29 +1935,77 @@ class YCone(Cone): """ _type = 'y-cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2') + + def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + R2 = kwargs.pop('R2', None) + if R2 is not None: + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) + r2 = R2 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) + + @property + def x0(self): + return self.coefficients['x0'] + + @property + def y0(self): + return self.coefficients['y0'] + + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r2(self): + return self.coefficients['r2'] + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r2.setter + def r2(self, r2): + check_type('r^2 coefficient', r2, Real) + self._coefficients['r2'] = r2 + + def _get_base_coeffs(self): + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + + b = -r2 + a = c = 1. + d = e = f = 0. + g, h, j = -2*x0, 2*y0*r2, -2*z0 + k = x0*x0 + z0*z0 - r2*y0*y0 + + return (a, b, c, d, e, f, g, h, j, k) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (z' - z_0)^2 - r^2(y' - y_0)^2` - - """ x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return x**2 + z**2 - self.r2*y**2 + return x*x + z*z - self.r2*y*y -class ZCone(Cone): +Cone.register(YCone) + + +class ZCone(QuadricMixin, Surface): """A cone parallel to the x-axis of the form :math:`(x - x_0)^2 + (y - y_0)^2 = r^2 (z - z_0)^2`. @@ -1877,29 +2054,77 @@ class ZCone(Cone): """ _type = 'z-cone' + _coeff_keys = ('x0', 'y0', 'z0', 'r2') + + def __init__(self, x0=0., y0=0., z0=0., r2=1., **kwargs): + R2 = kwargs.pop('R2', None) + if R2 is not None: + warn(_WARNING_UPPER.format(type(self).__name__, 'r2', 'R2'), + FutureWarning) + r2 = R2 + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (x0, y0, z0, r2)): + setattr(self, key, val) + + @property + def x0(self): + return self.coefficients['x0'] + + @property + def y0(self): + return self.coefficients['y0'] + + @property + def z0(self): + return self.coefficients['z0'] + + @property + def r2(self): + return self.coefficients['r2'] + + @x0.setter + def x0(self, x0): + check_type('x0 coefficient', x0, Real) + self._coefficients['x0'] = x0 + + @y0.setter + def y0(self, y0): + check_type('y0 coefficient', y0, Real) + self._coefficients['y0'] = y0 + + @z0.setter + def z0(self, z0): + check_type('z0 coefficient', z0, Real) + self._coefficients['z0'] = z0 + + @r2.setter + def r2(self, r2): + check_type('r^2 coefficient', r2, Real) + self._coefficients['r2'] = r2 + + def _get_base_coeffs(self): + x0, y0, z0, r2 = self.x0, self.y0, self.z0, self.r2 + + c = -r2 + a = b = 1. + d = e = f = 0. + g, h, j = -2*x0, -2*y0, 2*z0*r2 + k = x0*x0 + y0*y0 - r2*z0*z0 + + return (a, b, c, d, e, f, g, h, j, k) def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`(x' - x_0)^2 + (y' - y_0)^2 - r^2(z' - z_0)^2` - - """ x = point[0] - self.x0 y = point[1] - self.y0 z = point[2] - self.z0 - return x**2 + y**2 - self.r2*z**2 + return x*x + y*y - self.r2*z*z -class Quadric(Surface): +Cone.register(ZCone) + + +class Quadric(QuadricMixin, Surface): """A surface of the form :math:`Ax^2 + By^2 + Cz^2 + Dxy + Eyz + Fxz + Gx + Hy + Jz + K = 0`. @@ -1939,18 +2164,12 @@ class Quadric(Surface): _coeff_keys = ('a', 'b', 'c', 'd', 'e', 'f', 'g', 'h', 'j', 'k') def __init__(self, a=0., b=0., c=0., d=0., e=0., f=0., g=0., h=0., j=0., - k=0., boundary_type='transmission', name='', surface_id=None): - super().__init__(surface_id, boundary_type, name=name) - self.a = a - self.b = b - self.c = c - self.d = d - self.e = e - self.f = f - self.g = g - self.h = h - self.j = j - self.k = k + k=0., **kwargs): + + super().__init__(**kwargs) + + for key, val in zip(self._coeff_keys, (a, b, c, d, e, f, g, h, j, k)): + setattr(self, key, val) @property def a(self): @@ -2042,50 +2261,8 @@ class Quadric(Surface): check_type('k coefficient', k, Real) self._coefficients['k'] = k - def evaluate(self, point): - """Evaluate the surface equation at a given point. - - Parameters - ---------- - point : 3-tuple of float - The Cartesian coordinates, :math:`(x',y',z')`, at which the surface - equation should be evaluated. - - Returns - ------- - float - :math:`Ax'^2 + By'^2 + Cz'^2 + Dx'y' + Ey'z' + Fx'z' + Gx' + Hy' + - Jz' + K = 0` - - """ - x, y, z = point - return x*(self.a*x + self.d*y + self.g) + \ - y*(self.b*y + self.e*z + self.h) + \ - z*(self.c*z + self.f*x + self.j) + self.k - - def translate(self, vector): - """Translate surface in given direction - - Parameters - ---------- - vector : iterable of float - Direction in which surface should be translated - - Returns - ------- - openmc.Quadric - Translated surface - - """ - vx, vy, vz = vector - a, b, c, d, e, f, g, h, j, k = (getattr(self, key) for key in - self._coeff_keys) - k = (k + vx*vx + vy*vy + vz*vz + d*vx*vy + e*vy*vz + f*vx*vz - - g*vx - h*vy - j*vz) - g = g - 2*a*vx - d*vy - f*vz - h = h - 2*b*vy - d*vx - e*vz - j = j - 2*c*vz - e*vy - f*vx - return type(self)(a=a, b=b, c=c, d=d, e=e, f=f, g=g, h=h, j=j, k=k) + def _get_base_coeffs(self): + return tuple(getattr(self, c) for c in self._coeff_keys) class Halfspace(Region): @@ -2274,3 +2451,6 @@ class Halfspace(Region): # Return translated surface return type(self)(memo[key], self.side) + + +_SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()} diff --git a/src/cell.cpp b/src/cell.cpp index cd6199b755..1d175f05d2 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -377,7 +377,7 @@ CSGCell::CSGCell(pugi::xml_node cell_node) throw std::runtime_error{"Invalid surface ID " + std::to_string(abs(r)) + " specified in region for cell " + std::to_string(id_) + "."}; } - r = copysign(it->second + 1, r); + r = (r > 0) ? it->second + 1 : -(it->second + 1); } } @@ -536,8 +536,8 @@ CSGCell::to_hdf5(hid_t cell_group) const region_spec << " |"; } else { // Note the off-by-one indexing - region_spec << " " - << copysign(model::surfaces[abs(token)-1]->id_, token); + auto surf_id = model::surfaces[abs(token)-1]->id_; + region_spec << " " << ((token > 0) ? surf_id : -surf_id); } } write_string(group, "region", region_spec.str(), false); diff --git a/src/output.cpp b/src/output.cpp index d98cfd4a1d..18b94be9be 100644 --- a/src/output.cpp +++ b/src/output.cpp @@ -192,7 +192,7 @@ extern "C" void print_particle(Particle* p) // Display miscellaneous info. if (p->surface_ != 0) { const Surface& surf {*model::surfaces[std::abs(p->surface_)-1]}; - fmt::print(" Surface = {}\n", std::copysign(surf.id_, p->surface_)); + fmt::print(" Surface = {}\n", (p->surface_ > 0) ? surf.id_ : -surf.id_); } fmt::print(" Weight = {}\n", p->wgt_); if (settings::run_CE) { diff --git a/src/particle.cpp b/src/particle.cpp index d2d8b2e1b3..ea0d03cd01 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -1,7 +1,7 @@ #include "openmc/particle.h" #include // copy, min -#include // log, abs, copysign +#include // log, abs #include @@ -539,7 +539,7 @@ Particle::cross_surface() // TODO: off-by-one surface_ = rotational ? surf_p->i_periodic_ + 1 : - std::copysign(surf_p->i_periodic_ + 1, surface_); + ((surface_ > 0) ? surf_p->i_periodic_ + 1 : -(surf_p->i_periodic_ + 1)); // Figure out what cell particle is in now n_coord_ = 1; diff --git a/tests/unit_tests/test_data_photon.py b/tests/unit_tests/test_data_photon.py index c767d19e6c..f7274e9c19 100644 --- a/tests/unit_tests/test_data_photon.py +++ b/tests/unit_tests/test_data_photon.py @@ -2,6 +2,7 @@ from collections.abc import Mapping, Callable import os +from pathlib import Path import numpy as np import pandas as pd @@ -144,3 +145,9 @@ def test_export_to_hdf5(tmpdir, element): element2.bremsstrahlung['electron_energy']).all() # Export to hdf5 again element2.export_to_hdf5(filename, 'w') + +def test_photodat_only(run_in_tmpdir): + endf_dir = Path(os.environ['OPENMC_ENDF_DATA']) + photoatomic_file = endf_dir / 'photoat' / 'photoat-001_H_000.endf' + data = openmc.data.IncidentPhoton.from_endf(photoatomic_file) + data.export_to_hdf5('tmp.h5', 'w') \ No newline at end of file