mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 22:26:08 -04:00
Merge pull request #1491 from eepeterson/surface_refactoring2
Surface refactoring improvements
This commit is contained in:
Paul Romano
GitHub
commit
66dabca159
5 changed files with 397 additions and 169 deletions
|
|
@ -404,7 +404,7 @@ to install the Python package in :ref:`"editable" mode <devguide_editable>`.
|
|||
Prerequisites
|
||||
-------------
|
||||
|
||||
The Python API works with Python 3.4+. In addition to Python itself, the API
|
||||
The Python API works with Python 3.5+. In addition to Python itself, the API
|
||||
relies on a number of third-party packages. All prerequisites can be installed
|
||||
using Conda_ (recommended), pip_, or through the package manager in most Linux
|
||||
distributions. To run simulations in parallel using MPI, it is recommended to
|
||||
|
|
|
|||
|
|
@ -373,52 +373,7 @@ def get_hexagonal_prism(*args, **kwargs):
|
|||
return hexagonal_prism(*args, **kwargs)
|
||||
|
||||
|
||||
def cylinder_from_points(p1, p2, r, **kwargs):
|
||||
"""Return cylinder defined by two points passing through its center.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
p1, p2 : 3-tuples
|
||||
Coordinates of two points that pass through the center of the cylinder
|
||||
r : float
|
||||
Radius of the cylinder
|
||||
kwargs : dict
|
||||
Keyword arguments passed to the :class:`openmc.Quadric` constructor
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.Quadric
|
||||
Quadric surface representing the cylinder.
|
||||
|
||||
"""
|
||||
# Get x, y, z coordinates of two points
|
||||
x1, y1, z1 = p1
|
||||
x2, y2, z2 = p2
|
||||
|
||||
# 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.
|
||||
kwargs['a'] = dy*dy + dz*dz
|
||||
kwargs['b'] = dx*dx + dz*dz
|
||||
kwargs['c'] = dx*dx + dy*dy
|
||||
kwargs['d'] = -2*dx*dy
|
||||
kwargs['e'] = -2*dy*dz
|
||||
kwargs['f'] = -2*dx*dz
|
||||
kwargs['g'] = 2*(cy*dz - cz*dy)
|
||||
kwargs['h'] = 2*(cz*dx - cx*dz)
|
||||
kwargs['j'] = 2*(cx*dy - cy*dx)
|
||||
kwargs['k'] = cx*cx + cy*cy + cz*cz - (dx*dx + dy*dy + dz*dz)*r*r
|
||||
|
||||
return Quadric(**kwargs)
|
||||
cylinder_from_points = Cylinder.from_points
|
||||
|
||||
|
||||
def subdivide(surfaces):
|
||||
|
|
|
|||
|
|
@ -3,13 +3,14 @@ from collections import OrderedDict
|
|||
from copy import deepcopy
|
||||
from numbers import Real
|
||||
from xml.etree import ElementTree as ET
|
||||
from warnings import warn
|
||||
from warnings import warn, catch_warnings, simplefilter
|
||||
import math
|
||||
|
||||
import numpy as np
|
||||
|
||||
from openmc.checkvalue import check_type, check_value
|
||||
from openmc.region import Region, Intersection, Union
|
||||
from openmc.mixin import IDManagerMixin
|
||||
from openmc.mixin import IDManagerMixin, IDWarning
|
||||
|
||||
|
||||
_BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic', 'white']
|
||||
|
|
@ -361,6 +362,59 @@ class PlaneMixin(metaclass=ABCMeta):
|
|||
self._periodic_surface = periodic_surface
|
||||
periodic_surface._periodic_surface = self
|
||||
|
||||
def _get_base_coeffs(self):
|
||||
return (self.a, self.b, self.c, self.d)
|
||||
|
||||
def _get_normal(self):
|
||||
a, b, c = self._get_base_coeffs()[:3]
|
||||
return np.array((a, b, c)) / math.sqrt(a*a + b*b + c*c)
|
||||
|
||||
def bounding_box(self, side):
|
||||
"""Determine an axis-aligned bounding box.
|
||||
|
||||
An axis-aligned bounding box for Plane half-spaces is represented by
|
||||
its lower-left and upper-right coordinates. If the half-space is
|
||||
unbounded in a particular direction, numpy.inf is used to represent
|
||||
infinity.
|
||||
|
||||
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
|
||||
|
||||
"""
|
||||
# Compute the bounding box based on the normal vector to the plane
|
||||
nhat = self._get_normal()
|
||||
ll = np.array([-np.inf, -np.inf, -np.inf])
|
||||
ur = np.array([np.inf, np.inf, np.inf])
|
||||
# If the plane is axis aligned, find the proper bounding box
|
||||
if np.any(np.isclose(np.abs(nhat), 1., rtol=0., atol=self._atol)):
|
||||
sign = nhat.sum()
|
||||
a, b, c, d = self._get_base_coeffs()
|
||||
vals = [d/val if not np.isclose(val, 0., rtol=0., atol=self._atol)
|
||||
else np.nan for val in (a, b, c)]
|
||||
if side == '-':
|
||||
if sign > 0:
|
||||
ur = np.array([v if not np.isnan(v) else np.inf for v in vals])
|
||||
else:
|
||||
ll = np.array([v if not np.isnan(v) else -np.inf for v in vals])
|
||||
elif side == '+':
|
||||
if sign > 0:
|
||||
ll = np.array([v if not np.isnan(v) else -np.inf for v in vals])
|
||||
else:
|
||||
ur = np.array([v if not np.isnan(v) else np.inf for v in vals])
|
||||
|
||||
return (ll, ur)
|
||||
|
||||
def evaluate(self, point):
|
||||
"""Evaluate the surface equation at a given point.
|
||||
|
||||
|
|
@ -509,6 +563,12 @@ class Plane(PlaneMixin, Surface):
|
|||
FutureWarning)
|
||||
setattr(self, k.lower(), val)
|
||||
|
||||
@classmethod
|
||||
def __subclasshook__(cls, c):
|
||||
if cls is Plane and c in (XPlane, YPlane, ZPlane):
|
||||
return True
|
||||
return NotImplemented
|
||||
|
||||
@property
|
||||
def a(self):
|
||||
return self.coefficients['a']
|
||||
|
|
@ -545,9 +605,6 @@ class Plane(PlaneMixin, Surface):
|
|||
check_type('D coefficient', d, Real)
|
||||
self._coefficients['d'] = d
|
||||
|
||||
def _get_base_coeffs(self):
|
||||
return (self.a, self.b, self.c, self.d)
|
||||
|
||||
@classmethod
|
||||
def from_points(cls, p1, p2, p3, **kwargs):
|
||||
"""Return a plane given three points that pass through it.
|
||||
|
|
@ -638,29 +695,31 @@ class XPlane(PlaneMixin, Surface):
|
|||
def x0(self):
|
||||
return self.coefficients['x0']
|
||||
|
||||
@property
|
||||
def a(self):
|
||||
return 1.
|
||||
|
||||
@property
|
||||
def b(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def c(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def d(self):
|
||||
return self.x0
|
||||
|
||||
@x0.setter
|
||||
def x0(self, x0):
|
||||
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):
|
||||
if side == '-':
|
||||
return (np.array([-np.inf, -np.inf, -np.inf]),
|
||||
np.array([self.x0, np.inf, np.inf]))
|
||||
elif side == '+':
|
||||
return (np.array([self.x0, -np.inf, -np.inf]),
|
||||
np.array([np.inf, np.inf, np.inf]))
|
||||
|
||||
def evaluate(self, point):
|
||||
return point[0] - self.x0
|
||||
|
||||
|
||||
Plane.register(XPlane)
|
||||
|
||||
|
||||
class YPlane(PlaneMixin, Surface):
|
||||
"""A plane perpendicular to the y axis of the form :math:`y - y_0 = 0`
|
||||
|
||||
|
|
@ -718,29 +777,31 @@ class YPlane(PlaneMixin, Surface):
|
|||
def y0(self):
|
||||
return self.coefficients['y0']
|
||||
|
||||
@property
|
||||
def a(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def b(self):
|
||||
return 1.
|
||||
|
||||
@property
|
||||
def c(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def d(self):
|
||||
return self.y0
|
||||
|
||||
@y0.setter
|
||||
def y0(self, y0):
|
||||
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):
|
||||
if side == '-':
|
||||
return (np.array([-np.inf, -np.inf, -np.inf]),
|
||||
np.array([np.inf, self.y0, np.inf]))
|
||||
elif side == '+':
|
||||
return (np.array([-np.inf, self.y0, -np.inf]),
|
||||
np.array([np.inf, np.inf, np.inf]))
|
||||
|
||||
def evaluate(self, point):
|
||||
return point[1] - self.y0
|
||||
|
||||
|
||||
Plane.register(YPlane)
|
||||
|
||||
|
||||
class ZPlane(PlaneMixin, Surface):
|
||||
"""A plane perpendicular to the z axis of the form :math:`z - z_0 = 0`
|
||||
|
||||
|
|
@ -798,32 +859,43 @@ class ZPlane(PlaneMixin, Surface):
|
|||
def z0(self):
|
||||
return self.coefficients['z0']
|
||||
|
||||
@property
|
||||
def a(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def b(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def c(self):
|
||||
return 1.
|
||||
|
||||
@property
|
||||
def d(self):
|
||||
return self.z0
|
||||
|
||||
@z0.setter
|
||||
def z0(self, z0):
|
||||
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):
|
||||
if side == '-':
|
||||
return (np.array([-np.inf, -np.inf, -np.inf]),
|
||||
np.array([np.inf, np.inf, self.z0]))
|
||||
elif side == '+':
|
||||
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"""
|
||||
|
||||
@property
|
||||
def _origin(self):
|
||||
return np.array((self.x0, self.y0, self.z0))
|
||||
|
||||
@property
|
||||
def _axis(self):
|
||||
axis = np.array((self.dx, self.dy, self.dz))
|
||||
return axis / np.linalg.norm(axis)
|
||||
|
||||
def get_Abc(self, coeffs=None):
|
||||
"""Compute matrix, vector, and scalar coefficients for this surface or
|
||||
for a specified set of coefficients.
|
||||
|
|
@ -884,7 +956,7 @@ class QuadricMixin(metaclass=ABCMeta):
|
|||
"""
|
||||
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
|
||||
return x.T @ A @ x + b.T @ x + c
|
||||
|
||||
def translate(self, vector, inplace=False):
|
||||
"""Translate surface in given direction
|
||||
|
|
@ -907,17 +979,20 @@ class QuadricMixin(metaclass=ABCMeta):
|
|||
|
||||
surf = self if inplace else self.clone()
|
||||
|
||||
if set(('x0', 'y0', 'z0')).intersection(set(surf._coeff_keys)):
|
||||
if hasattr(self, 'x0'):
|
||||
for vi, xi in zip(vector, ('x0', 'y0', 'z0')):
|
||||
val = getattr(surf, xi, None)
|
||||
if val is not None:
|
||||
val = getattr(surf, xi)
|
||||
try:
|
||||
setattr(surf, xi, val + vi)
|
||||
except AttributeError:
|
||||
# That attribute is read only i.e x0 for XCylinder
|
||||
pass
|
||||
|
||||
else:
|
||||
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)
|
||||
g, h, j = bvec - 2*vector.T @ A
|
||||
k = cnst + vector.T @ A @ vector - bvec.T @ vector
|
||||
|
||||
for key, val in zip(('g', 'h', 'j', 'k'), (g, h, j, k)):
|
||||
setattr(surf, key, val)
|
||||
|
|
@ -993,16 +1068,17 @@ class Cylinder(QuadricMixin, Surface):
|
|||
_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)
|
||||
|
||||
@classmethod
|
||||
def __subclasshook__(cls, c):
|
||||
if cls is Cylinder and c in (XCylinder, YCylinder, ZCylinder):
|
||||
return True
|
||||
return NotImplemented
|
||||
|
||||
@property
|
||||
def x0(self):
|
||||
return self.coefficients['x0']
|
||||
|
|
@ -1068,8 +1144,8 @@ class Cylinder(QuadricMixin, Surface):
|
|||
|
||||
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
|
||||
x1, y1, z1 = self._origin
|
||||
x2, y2, z2 = self._origin + self._axis
|
||||
r = self.r
|
||||
|
||||
# Define intermediate terms
|
||||
|
|
@ -1118,10 +1194,6 @@ class Cylinder(QuadricMixin, Surface):
|
|||
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)
|
||||
|
|
@ -1130,6 +1202,24 @@ class Cylinder(QuadricMixin, Surface):
|
|||
|
||||
return cls(x0=x0, y0=y0, z0=z0, r=r, dx=dx, dy=dy, dz=dz, **kwargs)
|
||||
|
||||
def to_xml_element(self):
|
||||
"""Return XML representation of the surface
|
||||
|
||||
Returns
|
||||
-------
|
||||
element : xml.etree.ElementTree.Element
|
||||
XML element containing source data
|
||||
|
||||
"""
|
||||
# This method overrides Surface.to_xml_element to generate a Quadric
|
||||
# since the C++ layer doesn't support Cylinders right now
|
||||
with catch_warnings():
|
||||
simplefilter('ignore', IDWarning)
|
||||
kwargs = {'boundary_type': self.boundary_type, 'name': self.name,
|
||||
'surface_id': self.id}
|
||||
quad_rep = Quadric(*self._get_base_coeffs(), **kwargs)
|
||||
return quad_rep.to_xml_element()
|
||||
|
||||
|
||||
class XCylinder(QuadricMixin, Surface):
|
||||
"""An infinite cylinder whose length is parallel to the x-axis of the form
|
||||
|
|
@ -1202,6 +1292,22 @@ class XCylinder(QuadricMixin, Surface):
|
|||
def r(self):
|
||||
return self.coefficients['r']
|
||||
|
||||
@property
|
||||
def x0(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dx(self):
|
||||
return 1.
|
||||
|
||||
@property
|
||||
def dy(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dz(self):
|
||||
return 0.
|
||||
|
||||
@y0.setter
|
||||
def y0(self, y0):
|
||||
check_type('y0 coefficient', y0, Real)
|
||||
|
|
@ -1240,9 +1346,6 @@ class XCylinder(QuadricMixin, Surface):
|
|||
return y*y + z*z - self.r**2
|
||||
|
||||
|
||||
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`.
|
||||
|
|
@ -1314,6 +1417,22 @@ class YCylinder(QuadricMixin, Surface):
|
|||
def r(self):
|
||||
return self.coefficients['r']
|
||||
|
||||
@property
|
||||
def y0(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dx(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dy(self):
|
||||
return 1.
|
||||
|
||||
@property
|
||||
def dz(self):
|
||||
return 0.
|
||||
|
||||
@x0.setter
|
||||
def x0(self, x0):
|
||||
check_type('x0 coefficient', x0, Real)
|
||||
|
|
@ -1352,9 +1471,6 @@ class YCylinder(QuadricMixin, Surface):
|
|||
return x*x + z*z - self.r**2
|
||||
|
||||
|
||||
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`.
|
||||
|
|
@ -1426,6 +1542,22 @@ class ZCylinder(QuadricMixin, Surface):
|
|||
def r(self):
|
||||
return self.coefficients['r']
|
||||
|
||||
@property
|
||||
def z0(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dx(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dy(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dz(self):
|
||||
return 1.
|
||||
|
||||
@x0.setter
|
||||
def x0(self, x0):
|
||||
check_type('x0 coefficient', x0, Real)
|
||||
|
|
@ -1464,9 +1596,6 @@ class ZCylinder(QuadricMixin, Surface):
|
|||
return x*x + y*y - self.r**2
|
||||
|
||||
|
||||
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`.
|
||||
|
||||
|
|
@ -1656,9 +1785,6 @@ class Cone(QuadricMixin, Surface):
|
|||
_coeff_keys = ('x0', 'y0', 'z0', 'r2', 'dx', 'dy', 'dz')
|
||||
|
||||
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'),
|
||||
|
|
@ -1669,6 +1795,12 @@ class Cone(QuadricMixin, Surface):
|
|||
for key, val in zip(self._coeff_keys, (x0, y0, z0, r2, dx, dy, dz)):
|
||||
setattr(self, key, val)
|
||||
|
||||
@classmethod
|
||||
def __subclasshook__(cls, c):
|
||||
if cls is Cone and c in (XCone, YCone, ZCone):
|
||||
return True
|
||||
return NotImplemented
|
||||
|
||||
@property
|
||||
def x0(self):
|
||||
return self.coefficients['x0']
|
||||
|
|
@ -1745,27 +1877,42 @@ class Cone(QuadricMixin, Surface):
|
|||
# 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)
|
||||
x0, y0, z0 = self._origin
|
||||
dx, dy, dz = self._axis
|
||||
cos2 = 1 / (1 + self.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)
|
||||
a = cos2 - dx*dx
|
||||
b = cos2 - dy*dy
|
||||
c = cos2 - dz*dz
|
||||
d = -2*dx*dy
|
||||
e = -2*dy*dz
|
||||
f = -2*dx*dz
|
||||
g = 2*(dx*(dy*y0 + dz*z0) - a*x0)
|
||||
h = 2*(dy*(dx*x0 + dz*z0) - b*y0)
|
||||
j = 2*(dz*(dx*x0 + dy*y0) - c*z0)
|
||||
k = a*x0*x0 + b*y0*y0 + c*z0*z0 - 2*(dx*dy*x0*y0 + dy*dz*y0*z0 +
|
||||
dx*dz*x0*z0)
|
||||
|
||||
return (a, b, c, d, e, f, g, h, j, k)
|
||||
|
||||
def to_xml_element(self):
|
||||
"""Return XML representation of the surface
|
||||
|
||||
Returns
|
||||
-------
|
||||
element : xml.etree.ElementTree.Element
|
||||
XML element containing source data
|
||||
|
||||
"""
|
||||
# This method overrides Surface.to_xml_element to generate a Quadric
|
||||
# since the C++ layer doesn't support Cones right now
|
||||
with catch_warnings():
|
||||
simplefilter('ignore', IDWarning)
|
||||
kwargs = {'boundary_type': self.boundary_type, 'name': self.name,
|
||||
'surface_id': self.id}
|
||||
quad_rep = Quadric(*self._get_base_coeffs(), **kwargs)
|
||||
return quad_rep.to_xml_element()
|
||||
|
||||
|
||||
class XCone(QuadricMixin, Surface):
|
||||
"""A cone parallel to the x-axis of the form :math:`(y - y_0)^2 + (z - z_0)^2 =
|
||||
|
|
@ -1845,6 +1992,18 @@ class XCone(QuadricMixin, Surface):
|
|||
def r2(self):
|
||||
return self.coefficients['r2']
|
||||
|
||||
@property
|
||||
def dx(self):
|
||||
return 1.
|
||||
|
||||
@property
|
||||
def dy(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dz(self):
|
||||
return 0.
|
||||
|
||||
@x0.setter
|
||||
def x0(self, x0):
|
||||
check_type('x0 coefficient', x0, Real)
|
||||
|
|
@ -1883,9 +2042,6 @@ class XCone(QuadricMixin, Surface):
|
|||
return y*y + z*z - self.r2*x*x
|
||||
|
||||
|
||||
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`.
|
||||
|
|
@ -1964,6 +2120,18 @@ class YCone(QuadricMixin, Surface):
|
|||
def r2(self):
|
||||
return self.coefficients['r2']
|
||||
|
||||
@property
|
||||
def dx(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dy(self):
|
||||
return 1.
|
||||
|
||||
@property
|
||||
def dz(self):
|
||||
return 0.
|
||||
|
||||
@x0.setter
|
||||
def x0(self, x0):
|
||||
check_type('x0 coefficient', x0, Real)
|
||||
|
|
@ -2002,9 +2170,6 @@ class YCone(QuadricMixin, Surface):
|
|||
return x*x + z*z - self.r2*y*y
|
||||
|
||||
|
||||
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`.
|
||||
|
|
@ -2083,6 +2248,18 @@ class ZCone(QuadricMixin, Surface):
|
|||
def r2(self):
|
||||
return self.coefficients['r2']
|
||||
|
||||
@property
|
||||
def dx(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dy(self):
|
||||
return 0.
|
||||
|
||||
@property
|
||||
def dz(self):
|
||||
return 1.
|
||||
|
||||
@x0.setter
|
||||
def x0(self, x0):
|
||||
check_type('x0 coefficient', x0, Real)
|
||||
|
|
@ -2121,9 +2298,6 @@ class ZCone(QuadricMixin, Surface):
|
|||
return x*x + y*y - self.r2*z*z
|
||||
|
||||
|
||||
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`.
|
||||
|
|
@ -2452,5 +2626,4 @@ class Halfspace(Region):
|
|||
# Return translated surface
|
||||
return type(self)(memo[key], self.side)
|
||||
|
||||
|
||||
_SURFACE_CLASSES = {cls._type: cls for cls in Surface.__subclasses__()}
|
||||
|
|
|
|||
3
setup.py
3
setup.py
|
|
@ -56,14 +56,13 @@ kwargs = {
|
|||
'Topic :: Scientific/Engineering'
|
||||
'Programming Language :: C++',
|
||||
'Programming Language :: Python :: 3',
|
||||
'Programming Language :: Python :: 3.4',
|
||||
'Programming Language :: Python :: 3.5',
|
||||
'Programming Language :: Python :: 3.6',
|
||||
'Programming Language :: Python :: 3.7',
|
||||
],
|
||||
|
||||
# Dependencies
|
||||
'python_requires': '>=3.4',
|
||||
'python_requires': '>=3.5',
|
||||
'install_requires': [
|
||||
'numpy>=1.9', 'h5py', 'scipy', 'ipython', 'matplotlib',
|
||||
'pandas', 'lxml', 'uncertainties'
|
||||
|
|
|
|||
|
|
@ -2,6 +2,7 @@ from functools import partial
|
|||
from random import uniform, seed
|
||||
|
||||
import numpy as np
|
||||
import math
|
||||
import openmc
|
||||
import pytest
|
||||
|
||||
|
|
@ -139,6 +140,51 @@ def test_zplane():
|
|||
repr(s)
|
||||
|
||||
|
||||
def test_cylinder():
|
||||
x0, y0, z0, r = 2, 3, 4, 2
|
||||
dx, dy, dz = 1, -1, 1
|
||||
s = openmc.Cylinder(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r=r)
|
||||
assert s.x0 == 2
|
||||
assert s.y0 == 3
|
||||
assert s.z0 == 4
|
||||
assert s.dx == 1
|
||||
assert s.dy == -1
|
||||
assert s.dz == 1
|
||||
assert s.r == 2
|
||||
|
||||
# Check bounding box
|
||||
assert_infinite_bb(s)
|
||||
|
||||
# evaluate method
|
||||
# |(p - p1) ⨯ (p - p2)|^2 / |p2 - p1|^2 - r^2
|
||||
p1 = s._origin
|
||||
p2 = p1 + s._axis
|
||||
perp = np.array((1, -2, 1))*(1 / s._axis)
|
||||
divisor = np.linalg.norm(p2 - p1)
|
||||
pin = p1 + 5*s._axis # point inside cylinder
|
||||
pout = np.array((4., 0., 2.5)) # point outside the cylinder
|
||||
pon = p1 + s.r*perp / np.linalg.norm(perp) # point on cylinder
|
||||
for p, fn in zip((pin, pout, pon), (np.less, np.greater, np.isclose)):
|
||||
c1 = np.linalg.norm(np.cross(p - p1, p - p2)) / divisor
|
||||
val = c1*c1 - s.r*s.r
|
||||
p_eval = s.evaluate(p)
|
||||
assert fn(p_eval, 0.)
|
||||
assert p_eval == pytest.approx(val)
|
||||
|
||||
# translate method
|
||||
st = s.translate((1.0, 1.0, 1.0))
|
||||
assert st.x0 == s.x0 + 1
|
||||
assert st.y0 == s.y0 + 1
|
||||
assert st.z0 == s.z0 + 1
|
||||
assert st.dx == s.dx
|
||||
assert st.dy == s.dy
|
||||
assert st.dz == s.dz
|
||||
assert st.r == s.r
|
||||
|
||||
# Make sure repr works
|
||||
repr(s)
|
||||
|
||||
|
||||
def test_xcylinder():
|
||||
y, z, r = 3, 5, 2
|
||||
s = openmc.XCylinder(y0=y, z0=z, r=r)
|
||||
|
|
@ -284,6 +330,58 @@ def cone_common(apex, r2, cls):
|
|||
repr(s)
|
||||
|
||||
|
||||
def test_cone():
|
||||
x0, y0, z0, r2 = 2, 3, 4, 4
|
||||
dx, dy, dz = 1, -1, 1
|
||||
s = openmc.Cone(x0=x0, y0=y0, z0=z0, dx=dx, dy=dy, dz=dz, r2=r2)
|
||||
assert s.x0 == 2
|
||||
assert s.y0 == 3
|
||||
assert s.z0 == 4
|
||||
assert s.dx == 1
|
||||
assert s.dy == -1
|
||||
assert s.dz == 1
|
||||
assert s.r2 == 4
|
||||
|
||||
# Check bounding box
|
||||
assert_infinite_bb(s)
|
||||
|
||||
# evaluate method
|
||||
# cos^2(theta) * ((p - p1))**2 - (d @ (p - p1))^2
|
||||
# The argument r2 for cones is actually tan^2(theta) so that
|
||||
# cos^2(theta) = 1 / (1 + r2)
|
||||
#
|
||||
# This makes the evaluation equation shown below where p is the evaluation
|
||||
# point (x, y, z) p1 is the apex (origin) of the cone and r2 is related to
|
||||
# the aperature of the cone as described above
|
||||
# (p - p1) @ (p - p1) / (1 + r2) - (d @ (p - p1))^2
|
||||
# point inside
|
||||
p1 = s._origin
|
||||
d = s._axis
|
||||
perp = np.array((1, -2, 1))*(1 / d)
|
||||
perp /= np.linalg.norm(perp)
|
||||
pin = p1 + 5*d # point inside cone
|
||||
pout = p1 + 3.2*perp # point outside cone
|
||||
pon = p1 + 3.2*d + 3.2*math.sqrt(s.r2)*perp # point on cone
|
||||
for p, fn in zip((pin, pout, pon), (np.less, np.greater, np.isclose)):
|
||||
val = np.sum((p - p1)**2) / (1 + s.r2) - np.sum((d @ (p - p1))**2)
|
||||
p_eval = s.evaluate(p)
|
||||
assert fn(p_eval, 0.)
|
||||
assert p_eval == pytest.approx(val)
|
||||
|
||||
# translate method
|
||||
st = s.translate((1.0, 1.0, 1.0))
|
||||
assert st.x0 == s.x0 + 1
|
||||
assert st.y0 == s.y0 + 1
|
||||
assert st.z0 == s.z0 + 1
|
||||
assert st.dx == s.dx
|
||||
assert st.dy == s.dy
|
||||
assert st.dz == s.dz
|
||||
assert st.r2 == s.r2
|
||||
|
||||
# Make sure repr works
|
||||
repr(s)
|
||||
|
||||
|
||||
def test_xcone():
|
||||
apex = (10, 0, 0)
|
||||
r2 = 4
|
||||
|
|
@ -339,7 +437,7 @@ def test_cylinder_from_points():
|
|||
p1 = np.array([xi(), xi(), xi()])
|
||||
p2 = np.array([xi(), xi(), xi()])
|
||||
r = uniform(1.0, 100.0)
|
||||
s = openmc.model.cylinder_from_points(p1, p2, r)
|
||||
s = openmc.Cylinder.from_points(p1, p2, r)
|
||||
|
||||
# Points p1 and p2 need to be inside cylinder
|
||||
assert p1 in -s
|
||||
|
|
@ -369,24 +467,27 @@ def test_cylinder_from_points_axis():
|
|||
|
||||
# (x - 3)^2 + (y - 4)^2 = 2^2
|
||||
# x^2 + y^2 - 6x - 8y + 21 = 0
|
||||
s = openmc.model.cylinder_from_points((3., 4., 0.), (3., 4., 1.), 2.)
|
||||
assert (s.a, s.b, s.c) == pytest.approx((1., 1., 0.))
|
||||
assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
|
||||
assert (s.g, s.h, s.j) == pytest.approx((-6., -8., 0.))
|
||||
assert s.k == pytest.approx(21.)
|
||||
s = openmc.Cylinder.from_points((3., 4., 0.), (3., 4., 1.), 2.)
|
||||
a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs()
|
||||
assert (a, b, c) == pytest.approx((1., 1., 0.))
|
||||
assert (d, e, f) == pytest.approx((0., 0., 0.))
|
||||
assert (g, h, j) == pytest.approx((-6., -8., 0.))
|
||||
assert k == pytest.approx(21.)
|
||||
|
||||
# (y + 7)^2 + (z - 1)^2 = 3^2
|
||||
# y^2 + z^2 + 14y - 2z + 41 = 0
|
||||
s = openmc.model.cylinder_from_points((0., -7, 1.), (1., -7., 1.), 3.)
|
||||
assert (s.a, s.b, s.c) == pytest.approx((0., 1., 1.))
|
||||
assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
|
||||
assert (s.g, s.h, s.j) == pytest.approx((0., 14., -2.))
|
||||
assert s.k == 41.
|
||||
s = openmc.Cylinder.from_points((0., -7, 1.), (1., -7., 1.), 3.)
|
||||
a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs()
|
||||
assert (a, b, c) == pytest.approx((0., 1., 1.))
|
||||
assert (d, e, f) == pytest.approx((0., 0., 0.))
|
||||
assert (g, h, j) == pytest.approx((0., 14., -2.))
|
||||
assert k == 41.
|
||||
|
||||
# (x - 2)^2 + (z - 5)^2 = 4^2
|
||||
# x^2 + z^2 - 4x - 10z + 13 = 0
|
||||
s = openmc.model.cylinder_from_points((2., 0., 5.), (2., 1., 5.), 4.)
|
||||
assert (s.a, s.b, s.c) == pytest.approx((1., 0., 1.))
|
||||
assert (s.d, s.e, s.f) == pytest.approx((0., 0., 0.))
|
||||
assert (s.g, s.h, s.j) == pytest.approx((-4., 0., -10.))
|
||||
assert s.k == pytest.approx(13.)
|
||||
s = openmc.Cylinder.from_points((2., 0., 5.), (2., 1., 5.), 4.)
|
||||
a, b, c, d, e, f, g, h, j, k = s._get_base_coeffs()
|
||||
assert (a, b, c) == pytest.approx((1., 0., 1.))
|
||||
assert (d, e, f) == pytest.approx((0., 0., 0.))
|
||||
assert (g, h, j) == pytest.approx((-4., 0., -10.))
|
||||
assert k == pytest.approx(13.)
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue