mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 06:05:58 -04:00
Move get_*_prism functions to openmc.model
This commit is contained in:
Paul Romano
parent
198df07311
commit
6e9b8a5ce5
5 changed files with 266 additions and 263 deletions
|
|
@ -91,18 +91,6 @@ Many of the above classes are derived from several abstract classes:
|
|||
openmc.Region
|
||||
openmc.Lattice
|
||||
|
||||
Two helper function are also available to create rectangular and hexagonal
|
||||
prisms defined by the intersection of four and six surface half-spaces,
|
||||
respectively.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated
|
||||
:nosignatures:
|
||||
:template: myfunction.rst
|
||||
|
||||
openmc.get_hexagonal_prism
|
||||
openmc.get_rectangular_prism
|
||||
|
||||
.. _pythonapi_tallies:
|
||||
|
||||
Constructing Tallies
|
||||
|
|
|
|||
|
|
@ -5,11 +5,18 @@
|
|||
Convenience Functions
|
||||
---------------------
|
||||
|
||||
Several helper functions are available here. Ther first two create rectangular
|
||||
and hexagonal prisms defined by the intersection of four and six surface
|
||||
half-spaces, respectively. The last function takes a sequence of surfaces and
|
||||
returns the regions that separate them.
|
||||
|
||||
.. autosummary::
|
||||
:toctree: generated
|
||||
:nosignatures:
|
||||
:template: myfunction.rst
|
||||
|
||||
openmc.model.get_hexagonal_prism
|
||||
openmc.model.get_rectangular_prism
|
||||
openmc.model.subdivide
|
||||
|
||||
TRISO Fuel Modeling
|
||||
|
|
|
|||
|
|
@ -28,4 +28,7 @@ from openmc.mixin import *
|
|||
from openmc.plotter import *
|
||||
from openmc.search import *
|
||||
|
||||
# Import a few convencience functions that used to be here
|
||||
from openmc.model import get_rectangular_prism, get_hexagonal_prism
|
||||
|
||||
__version__ = '0.9.0'
|
||||
|
|
|
|||
|
|
@ -1,3 +1,257 @@
|
|||
from __future__ import division
|
||||
from collections import Iterable, OrderedDict
|
||||
from math import sqrt
|
||||
from numbers import Real
|
||||
|
||||
from openmc import XPlane, YPlane, Plane, ZCylinder
|
||||
from openmc.checkvalue import check_type, check_value
|
||||
|
||||
|
||||
def get_rectangular_prism(width, height, axis='z', origin=(0., 0.),
|
||||
boundary_type='transmission', corner_radius=0.):
|
||||
"""Get an infinite rectangular prism from four planar surfaces.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
width: float
|
||||
Prism width in units of cm. The width is aligned with the y, x,
|
||||
or x axes for prisms parallel to the x, y, or z axis, respectively.
|
||||
height: float
|
||||
Prism height in units of cm. The height is aligned with the z, z,
|
||||
or y axes for prisms parallel to the x, y, or z axis, respectively.
|
||||
axis : {'x', 'y', 'z'}
|
||||
Axis with which the infinite length of the prism should be aligned.
|
||||
Defaults to 'z'.
|
||||
origin: Iterable of two floats
|
||||
Origin of the prism. The two floats correspond to (y,z), (x,z) or
|
||||
(x,y) for prisms parallel to the x, y or z axis, respectively.
|
||||
Defaults to (0., 0.).
|
||||
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
|
||||
Boundary condition that defines the behavior for particles hitting the
|
||||
surfaces comprising the rectangular prism (default is 'transmission').
|
||||
corner_radius: float
|
||||
Prism corner radius in units of cm. Defaults to 0.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.Region
|
||||
The inside of a rectangular prism
|
||||
|
||||
"""
|
||||
|
||||
check_type('width', width, Real)
|
||||
check_type('height', height, Real)
|
||||
check_type('corner_radius', corner_radius, Real)
|
||||
check_value('axis', axis, ['x', 'y', 'z'])
|
||||
check_type('origin', origin, Iterable, Real)
|
||||
|
||||
# Define function to create a plane on given axis
|
||||
def plane(axis, name, value):
|
||||
cls = globals()['{}Plane'.format(axis.upper())]
|
||||
return cls(name='{} {}'.format(name, axis),
|
||||
boundary_type=boundary_type,
|
||||
**{axis + '0': value})
|
||||
|
||||
if axis == 'x':
|
||||
x1, x2 = 'y', 'z'
|
||||
elif axis == 'y':
|
||||
x1, x2 = 'x', 'z'
|
||||
else:
|
||||
x1, x2 = 'x', 'y'
|
||||
|
||||
# Get cylinder class corresponding to given axis
|
||||
cyl = globals()['{}Cylinder'.format(axis.upper())]
|
||||
|
||||
# Create rectangular region
|
||||
min_x1 = plane(x1, 'minimum', -width/2 + origin[0])
|
||||
max_x1 = plane(x1, 'maximum', width/2 + origin[0])
|
||||
min_x2 = plane(x2, 'minimum', -height/2 + origin[1])
|
||||
max_x2 = plane(x2, 'maximum', height/2 + origin[1])
|
||||
prism = +min_x1 & -max_x1 & +min_x2 & -max_x2
|
||||
|
||||
# Handle rounded corners if given
|
||||
if corner_radius > 0.:
|
||||
args = {'R': corner_radius, 'boundary_type': boundary_type}
|
||||
|
||||
args[x1 + '0'] = origin[0] - width/2 + corner_radius
|
||||
args[x2 + '0'] = origin[1] - height/2 + corner_radius
|
||||
x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args)
|
||||
|
||||
args[x1 + '0'] = origin[0] - width/2 + corner_radius
|
||||
args[x2 + '0'] = origin[1] - height/2 + corner_radius
|
||||
x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args)
|
||||
|
||||
args[x1 + '0'] = origin[0] - width/2 + corner_radius
|
||||
args[x2 + '0'] = origin[1] + height/2 - corner_radius
|
||||
x1_min_x2_max = cyl(name='{} min {} max'.format(x1, x2), **args)
|
||||
|
||||
args[x1 + '0'] = origin[0] + width/2 - corner_radius
|
||||
args[x2 + '0'] = origin[1] - height/2 + corner_radius
|
||||
x1_max_x2_min = cyl(name='{} max {} min'.format(x1, x2), **args)
|
||||
|
||||
args[x1 + '0'] = origin[0] + width/2 - corner_radius
|
||||
args[x2 + '0'] = origin[1] + height/2 - corner_radius
|
||||
x1_max_x2_max = cyl(name='{} max {} max'.format(x1, x2), **args)
|
||||
|
||||
x1_min = plane(x1, 'min', -width/2 + origin[0] + corner_radius)
|
||||
x1_max = plane(x1, 'max', width/2 + origin[0] - corner_radius)
|
||||
x2_min = plane(x2, 'min', -height/2 + origin[1] + corner_radius)
|
||||
x2_max = plane(x2, 'max', height/2 + origin[1] - corner_radius)
|
||||
|
||||
corners = (+x1_min_x2_min & -x1_min & -x2_min) | \
|
||||
(+x1_min_x2_max & -x1_min & +x2_max) | \
|
||||
(+x1_max_x2_min & +x1_max & -x2_min) | \
|
||||
(+x1_max_x2_max & +x1_max & +x2_max)
|
||||
|
||||
prism = prism & ~corners
|
||||
|
||||
return prism
|
||||
|
||||
|
||||
def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
|
||||
boundary_type='transmission', corner_radius=0.):
|
||||
"""Create a hexagon region from six surface planes.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
edge_length : float
|
||||
Length of a side of the hexagon in cm
|
||||
orientation : {'x', 'y'}
|
||||
An 'x' orientation means that two sides of the hexagon are parallel to
|
||||
the x-axis and a 'y' orientation means that two sides of the hexagon are
|
||||
parallel to the y-axis.
|
||||
origin: Iterable of two floats
|
||||
Origin of the prism. Defaults to (0., 0.).
|
||||
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
|
||||
Boundary condition that defines the behavior for particles hitting the
|
||||
surfaces comprising the hexagonal prism (default is 'transmission').
|
||||
corner_radius: float
|
||||
Prism corner radius in units of cm. Defaults to 0.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.Region
|
||||
The inside of a hexagonal prism
|
||||
|
||||
"""
|
||||
|
||||
l = edge_length
|
||||
x, y = origin
|
||||
|
||||
if orientation == 'y':
|
||||
right = XPlane(x0=x + sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
left = XPlane(x0=x - sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
c = sqrt(3.)/3.
|
||||
|
||||
# y = -x/sqrt(3) + a
|
||||
upper_right = Plane(A=c, B=1., D=l+x*c+y, boundary_type=boundary_type)
|
||||
|
||||
# y = x/sqrt(3) + a
|
||||
upper_left = Plane(A=-c, B=1., D=l-x*c+y, boundary_type=boundary_type)
|
||||
|
||||
# y = x/sqrt(3) - a
|
||||
lower_right = Plane(A=-c, B=1., D=-l-x*c+y, boundary_type=boundary_type)
|
||||
|
||||
# y = -x/sqrt(3) - a
|
||||
lower_left = Plane(A=c, B=1., D=-l+x*c+y, boundary_type=boundary_type)
|
||||
|
||||
prism = -right & +left & -upper_right & -upper_left & \
|
||||
+lower_right & +lower_left
|
||||
|
||||
if boundary_type == 'periodic':
|
||||
right.periodic_surface = left
|
||||
upper_right.periodic_surface = lower_left
|
||||
lower_right.periodic_surface = upper_left
|
||||
|
||||
elif orientation == 'x':
|
||||
top = YPlane(y0=y + sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
bottom = YPlane(y0=y - sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
c = sqrt(3.)
|
||||
|
||||
# y = -sqrt(3)*(x - a)
|
||||
upper_right = Plane(A=c, B=1., D=c*l+x*c+y, boundary_type=boundary_type)
|
||||
|
||||
# y = sqrt(3)*(x + a)
|
||||
lower_right = Plane(A=-c, B=1., D=-c*l-x*c+y,
|
||||
boundary_type=boundary_type)
|
||||
|
||||
# y = -sqrt(3)*(x + a)
|
||||
lower_left = Plane(A=c, B=1., D=-c*l+x*c+y, boundary_type=boundary_type)
|
||||
|
||||
# y = sqrt(3)*(x + a)
|
||||
upper_left = Plane(A=-c, B=1., D=c*l-x*c+y, boundary_type=boundary_type)
|
||||
|
||||
prism = -top & +bottom & -upper_right & +lower_right & \
|
||||
+lower_left & -upper_left
|
||||
|
||||
if boundary_type == 'periodic':
|
||||
top.periodic_surface = bottom
|
||||
upper_right.periodic_surface = lower_left
|
||||
lower_right.periodic_surface = upper_left
|
||||
|
||||
# Handle rounded corners if given
|
||||
if corner_radius > 0.:
|
||||
if boundary_type == 'periodic':
|
||||
raise ValueError('Periodic boundary conditions not permitted when '
|
||||
'rounded corners are used.')
|
||||
|
||||
c = sqrt(3.)/2
|
||||
t = l - corner_radius/c
|
||||
|
||||
# Cylinder with corner radius and boundary type pre-applied
|
||||
cyl1 = partial(ZCylinder, R=corner_radius, boundary_type=boundary_type)
|
||||
cyl2 = partial(ZCylinder, R=corner_radius/(2*c),
|
||||
boundary_type=boundary_type)
|
||||
|
||||
if orientation == 'x':
|
||||
x_min_y_min_in = cyl1(name='x min y min in', x0=x-t/2, y0=y-c*t)
|
||||
x_min_y_max_in = cyl1(name='x min y max in', x0=x+t/2, y0=y-c*t)
|
||||
x_max_y_min_in = cyl1(name='x max y min in', x0=x-t/2, y0=y+c*t)
|
||||
x_max_y_max_in = cyl1(name='x max y max in', x0=x+t/2, y0=y+c*t)
|
||||
x_min_in = cyl1(name='x min in', x0=x-t, y0=y)
|
||||
x_max_in = cyl1(name='x max in', x0=x+t, y0=y)
|
||||
|
||||
x_min_y_min_out = cyl2(name='x min y min out', x0=x-l/2, y0=y-c*l)
|
||||
x_min_y_max_out = cyl2(name='x min y max out', x0=x+l/2, y0=y-c*l)
|
||||
x_max_y_min_out = cyl2(name='x max y min out', x0=x-l/2, y0=y+c*l)
|
||||
x_max_y_max_out = cyl2(name='x max y max out', x0=x+l/2, y0=y+c*l)
|
||||
x_min_out = cyl2(name='x min out', x0=x-l, y0=y)
|
||||
x_max_out = cyl2(name='x max out', x0=x+l, y0=y)
|
||||
|
||||
corners = (+x_min_y_min_in & -x_min_y_min_out |
|
||||
+x_min_y_max_in & -x_min_y_max_out |
|
||||
+x_max_y_min_in & -x_max_y_min_out |
|
||||
+x_max_y_max_in & -x_max_y_max_out |
|
||||
+x_min_in & -x_min_out |
|
||||
+x_max_in & -x_max_out)
|
||||
|
||||
elif orientation == 'y':
|
||||
x_min_y_min_in = cyl1(name='x min y min in', x0=x-c*t, y0=y-t/2)
|
||||
x_min_y_max_in = cyl1(name='x min y max in', x0=x-c*t, y0=y+t/2)
|
||||
x_max_y_min_in = cyl1(name='x max y min in', x0=x+c*t, y0=y-t/2)
|
||||
x_max_y_max_in = cyl1(name='x max y max in', x0=x+c*t, y0=y+t/2)
|
||||
y_min_in = cyl1(name='y min in', x0=x, y0=y-t)
|
||||
y_max_in = cyl1(name='y max in', x0=x, y0=y+t)
|
||||
|
||||
x_min_y_min_out = cyl2(name='x min y min out', x0=x-c*l, y0=y-l/2)
|
||||
x_min_y_max_out = cyl2(name='x min y max out', x0=x-c*l, y0=y+l/2)
|
||||
x_max_y_min_out = cyl2(name='x max y min out', x0=x+c*l, y0=y-l/2)
|
||||
x_max_y_max_out = cyl2(name='x max y max out', x0=x+c*l, y0=y+l/2)
|
||||
y_min_out = cyl2(name='y min out', x0=x, y0=y-l)
|
||||
y_max_out = cyl2(name='y max out', x0=x, y0=y+l)
|
||||
|
||||
corners = (+x_min_y_min_in & -x_min_y_min_out |
|
||||
+x_min_y_max_in & -x_min_y_max_out |
|
||||
+x_max_y_min_in & -x_max_y_min_out |
|
||||
+x_max_y_max_in & -x_max_y_max_out |
|
||||
+y_min_in & -y_min_out |
|
||||
+y_max_in & -y_max_out)
|
||||
|
||||
prism = prism & ~corners
|
||||
|
||||
return prism
|
||||
|
||||
|
||||
def subdivide(surfaces):
|
||||
"""Create regions separated by a series of surfaces.
|
||||
|
||||
|
|
|
|||
|
|
@ -1,23 +1,19 @@
|
|||
from __future__ import division
|
||||
from abc import ABCMeta
|
||||
from collections import Iterable, OrderedDict
|
||||
from collections import OrderedDict
|
||||
from copy import deepcopy
|
||||
from functools import partial
|
||||
from numbers import Real, Integral
|
||||
from xml.etree import ElementTree as ET
|
||||
from math import sqrt
|
||||
|
||||
from six import add_metaclass, string_types
|
||||
import numpy as np
|
||||
|
||||
from openmc.checkvalue import check_type, check_value, check_greater_than
|
||||
from openmc.checkvalue import check_type, check_value
|
||||
from openmc.region import Region, Intersection, Union
|
||||
from openmc.mixin import IDManagerMixin
|
||||
|
||||
|
||||
# A static variable for auto-generated Surface IDs
|
||||
AUTO_SURFACE_ID = 10000
|
||||
|
||||
_BOUNDARY_TYPES = ['transmission', 'vacuum', 'reflective', 'periodic']
|
||||
|
||||
|
||||
|
|
@ -1944,248 +1940,3 @@ class Halfspace(Region):
|
|||
clone = deepcopy(self)
|
||||
clone.surface = self.surface.clone(memo)
|
||||
return clone
|
||||
|
||||
|
||||
def get_rectangular_prism(width, height, axis='z', origin=(0., 0.),
|
||||
boundary_type='transmission', corner_radius=0.):
|
||||
"""Get an infinite rectangular prism from four planar surfaces.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
width: float
|
||||
Prism width in units of cm. The width is aligned with the y, x,
|
||||
or x axes for prisms parallel to the x, y, or z axis, respectively.
|
||||
height: float
|
||||
Prism height in units of cm. The height is aligned with the z, z,
|
||||
or y axes for prisms parallel to the x, y, or z axis, respectively.
|
||||
axis : {'x', 'y', 'z'}
|
||||
Axis with which the infinite length of the prism should be aligned.
|
||||
Defaults to 'z'.
|
||||
origin: Iterable of two floats
|
||||
Origin of the prism. The two floats correspond to (y,z), (x,z) or
|
||||
(x,y) for prisms parallel to the x, y or z axis, respectively.
|
||||
Defaults to (0., 0.).
|
||||
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
|
||||
Boundary condition that defines the behavior for particles hitting the
|
||||
surfaces comprising the rectangular prism (default is 'transmission').
|
||||
corner_radius: float
|
||||
Prism corner radius in units of cm. Defaults to 0.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.Region
|
||||
The inside of a rectangular prism
|
||||
|
||||
"""
|
||||
|
||||
check_type('width', width, Real)
|
||||
check_type('height', height, Real)
|
||||
check_type('corner_radius', corner_radius, Real)
|
||||
check_value('axis', axis, ['x', 'y', 'z'])
|
||||
check_type('origin', origin, Iterable, Real)
|
||||
|
||||
# Define function to create a plane on given axis
|
||||
def plane(axis, name, value):
|
||||
cls = globals()['{}Plane'.format(axis.upper())]
|
||||
return cls(name='{} {}'.format(name, axis),
|
||||
boundary_type=boundary_type,
|
||||
**{axis + '0': value})
|
||||
|
||||
if axis == 'x':
|
||||
x1, x2 = 'y', 'z'
|
||||
elif axis == 'y':
|
||||
x1, x2 = 'x', 'z'
|
||||
else:
|
||||
x1, x2 = 'x', 'y'
|
||||
|
||||
# Get cylinder class corresponding to given axis
|
||||
cyl = globals()['{}Cylinder'.format(axis.upper())]
|
||||
|
||||
# Create rectangular region
|
||||
min_x1 = plane(x1, 'minimum', -width/2 + origin[0])
|
||||
max_x1 = plane(x1, 'maximum', width/2 + origin[0])
|
||||
min_x2 = plane(x2, 'minimum', -height/2 + origin[1])
|
||||
max_x2 = plane(x2, 'maximum', height/2 + origin[1])
|
||||
prism = +min_x1 & -max_x1 & +min_x2 & -max_x2
|
||||
|
||||
# Handle rounded corners if given
|
||||
if corner_radius > 0.:
|
||||
args = {'R': corner_radius, 'boundary_type': boundary_type}
|
||||
|
||||
args[x1 + '0'] = origin[0] - width/2 + corner_radius
|
||||
args[x2 + '0'] = origin[1] - height/2 + corner_radius
|
||||
x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args)
|
||||
|
||||
args[x1 + '0'] = origin[0] - width/2 + corner_radius
|
||||
args[x2 + '0'] = origin[1] - height/2 + corner_radius
|
||||
x1_min_x2_min = cyl(name='{} min {} min'.format(x1, x2), **args)
|
||||
|
||||
args[x1 + '0'] = origin[0] - width/2 + corner_radius
|
||||
args[x2 + '0'] = origin[1] + height/2 - corner_radius
|
||||
x1_min_x2_max = cyl(name='{} min {} max'.format(x1, x2), **args)
|
||||
|
||||
args[x1 + '0'] = origin[0] + width/2 - corner_radius
|
||||
args[x2 + '0'] = origin[1] - height/2 + corner_radius
|
||||
x1_max_x2_min = cyl(name='{} max {} min'.format(x1, x2), **args)
|
||||
|
||||
args[x1 + '0'] = origin[0] + width/2 - corner_radius
|
||||
args[x2 + '0'] = origin[1] + height/2 - corner_radius
|
||||
x1_max_x2_max = cyl(name='{} max {} max'.format(x1, x2), **args)
|
||||
|
||||
x1_min = plane(x1, 'min', -width/2 + origin[0] + corner_radius)
|
||||
x1_max = plane(x1, 'max', width/2 + origin[0] - corner_radius)
|
||||
x2_min = plane(x2, 'min', -height/2 + origin[1] + corner_radius)
|
||||
x2_max = plane(x2, 'max', height/2 + origin[1] - corner_radius)
|
||||
|
||||
corners = (+x1_min_x2_min & -x1_min & -x2_min) | \
|
||||
(+x1_min_x2_max & -x1_min & +x2_max) | \
|
||||
(+x1_max_x2_min & +x1_max & -x2_min) | \
|
||||
(+x1_max_x2_max & +x1_max & +x2_max)
|
||||
|
||||
prism = prism & ~corners
|
||||
|
||||
return prism
|
||||
|
||||
|
||||
def get_hexagonal_prism(edge_length=1., orientation='y', origin=(0., 0.),
|
||||
boundary_type='transmission', corner_radius=0.):
|
||||
"""Create a hexagon region from six surface planes.
|
||||
|
||||
Parameters
|
||||
----------
|
||||
edge_length : float
|
||||
Length of a side of the hexagon in cm
|
||||
orientation : {'x', 'y'}
|
||||
An 'x' orientation means that two sides of the hexagon are parallel to
|
||||
the x-axis and a 'y' orientation means that two sides of the hexagon are
|
||||
parallel to the y-axis.
|
||||
origin: Iterable of two floats
|
||||
Origin of the prism. Defaults to (0., 0.).
|
||||
boundary_type : {'transmission, 'vacuum', 'reflective', 'periodic'}
|
||||
Boundary condition that defines the behavior for particles hitting the
|
||||
surfaces comprising the hexagonal prism (default is 'transmission').
|
||||
corner_radius: float
|
||||
Prism corner radius in units of cm. Defaults to 0.
|
||||
|
||||
Returns
|
||||
-------
|
||||
openmc.Region
|
||||
The inside of a hexagonal prism
|
||||
|
||||
"""
|
||||
|
||||
l = edge_length
|
||||
x, y = origin
|
||||
|
||||
if orientation == 'y':
|
||||
right = XPlane(x0=x + sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
left = XPlane(x0=x - sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
c = sqrt(3.)/3.
|
||||
|
||||
# y = -x/sqrt(3) + a
|
||||
upper_right = Plane(A=c, B=1., D=l+x*c+y, boundary_type=boundary_type)
|
||||
|
||||
# y = x/sqrt(3) + a
|
||||
upper_left = Plane(A=-c, B=1., D=l-x*c+y, boundary_type=boundary_type)
|
||||
|
||||
# y = x/sqrt(3) - a
|
||||
lower_right = Plane(A=-c, B=1., D=-l-x*c+y, boundary_type=boundary_type)
|
||||
|
||||
# y = -x/sqrt(3) - a
|
||||
lower_left = Plane(A=c, B=1., D=-l+x*c+y, boundary_type=boundary_type)
|
||||
|
||||
prism = -right & +left & -upper_right & -upper_left & \
|
||||
+lower_right & +lower_left
|
||||
|
||||
if boundary_type == 'periodic':
|
||||
right.periodic_surface = left
|
||||
upper_right.periodic_surface = lower_left
|
||||
lower_right.periodic_surface = upper_left
|
||||
|
||||
elif orientation == 'x':
|
||||
top = YPlane(y0=y + sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
bottom = YPlane(y0=y - sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
c = sqrt(3.)
|
||||
|
||||
# y = -sqrt(3)*(x - a)
|
||||
upper_right = Plane(A=c, B=1., D=c*l+x*c+y, boundary_type=boundary_type)
|
||||
|
||||
# y = sqrt(3)*(x + a)
|
||||
lower_right = Plane(A=-c, B=1., D=-c*l-x*c+y,
|
||||
boundary_type=boundary_type)
|
||||
|
||||
# y = -sqrt(3)*(x + a)
|
||||
lower_left = Plane(A=c, B=1., D=-c*l+x*c+y, boundary_type=boundary_type)
|
||||
|
||||
# y = sqrt(3)*(x + a)
|
||||
upper_left = Plane(A=-c, B=1., D=c*l-x*c+y, boundary_type=boundary_type)
|
||||
|
||||
prism = -top & +bottom & -upper_right & +lower_right & \
|
||||
+lower_left & -upper_left
|
||||
|
||||
if boundary_type == 'periodic':
|
||||
top.periodic_surface = bottom
|
||||
upper_right.periodic_surface = lower_left
|
||||
lower_right.periodic_surface = upper_left
|
||||
|
||||
# Handle rounded corners if given
|
||||
if corner_radius > 0.:
|
||||
if boundary_type == 'periodic':
|
||||
raise ValueError('Periodic boundary conditions not permitted when '
|
||||
'rounded corners are used.')
|
||||
|
||||
c = sqrt(3.)/2
|
||||
t = l - corner_radius/c
|
||||
|
||||
# Cylinder with corner radius and boundary type pre-applied
|
||||
cyl1 = partial(ZCylinder, R=corner_radius, boundary_type=boundary_type)
|
||||
cyl2 = partial(ZCylinder, R=corner_radius/(2*c),
|
||||
boundary_type=boundary_type)
|
||||
|
||||
if orientation == 'x':
|
||||
x_min_y_min_in = cyl1(name='x min y min in', x0=x-t/2, y0=y-c*t)
|
||||
x_min_y_max_in = cyl1(name='x min y max in', x0=x+t/2, y0=y-c*t)
|
||||
x_max_y_min_in = cyl1(name='x max y min in', x0=x-t/2, y0=y+c*t)
|
||||
x_max_y_max_in = cyl1(name='x max y max in', x0=x+t/2, y0=y+c*t)
|
||||
x_min_in = cyl1(name='x min in', x0=x-t, y0=y)
|
||||
x_max_in = cyl1(name='x max in', x0=x+t, y0=y)
|
||||
|
||||
x_min_y_min_out = cyl2(name='x min y min out', x0=x-l/2, y0=y-c*l)
|
||||
x_min_y_max_out = cyl2(name='x min y max out', x0=x+l/2, y0=y-c*l)
|
||||
x_max_y_min_out = cyl2(name='x max y min out', x0=x-l/2, y0=y+c*l)
|
||||
x_max_y_max_out = cyl2(name='x max y max out', x0=x+l/2, y0=y+c*l)
|
||||
x_min_out = cyl2(name='x min out', x0=x-l, y0=y)
|
||||
x_max_out = cyl2(name='x max out', x0=x+l, y0=y)
|
||||
|
||||
corners = (+x_min_y_min_in & -x_min_y_min_out |
|
||||
+x_min_y_max_in & -x_min_y_max_out |
|
||||
+x_max_y_min_in & -x_max_y_min_out |
|
||||
+x_max_y_max_in & -x_max_y_max_out |
|
||||
+x_min_in & -x_min_out |
|
||||
+x_max_in & -x_max_out)
|
||||
|
||||
elif orientation == 'y':
|
||||
x_min_y_min_in = cyl1(name='x min y min in', x0=x-c*t, y0=y-t/2)
|
||||
x_min_y_max_in = cyl1(name='x min y max in', x0=x-c*t, y0=y+t/2)
|
||||
x_max_y_min_in = cyl1(name='x max y min in', x0=x+c*t, y0=y-t/2)
|
||||
x_max_y_max_in = cyl1(name='x max y max in', x0=x+c*t, y0=y+t/2)
|
||||
y_min_in = cyl1(name='y min in', x0=x, y0=y-t)
|
||||
y_max_in = cyl1(name='y max in', x0=x, y0=y+t)
|
||||
|
||||
x_min_y_min_out = cyl2(name='x min y min out', x0=x-c*l, y0=y-l/2)
|
||||
x_min_y_max_out = cyl2(name='x min y max out', x0=x-c*l, y0=y+l/2)
|
||||
x_max_y_min_out = cyl2(name='x max y min out', x0=x+c*l, y0=y-l/2)
|
||||
x_max_y_max_out = cyl2(name='x max y max out', x0=x+c*l, y0=y+l/2)
|
||||
y_min_out = cyl2(name='y min out', x0=x, y0=y-l)
|
||||
y_max_out = cyl2(name='y max out', x0=x, y0=y+l)
|
||||
|
||||
corners = (+x_min_y_min_in & -x_min_y_min_out |
|
||||
+x_min_y_max_in & -x_min_y_max_out |
|
||||
+x_max_y_min_in & -x_max_y_min_out |
|
||||
+x_max_y_max_in & -x_max_y_max_out |
|
||||
+y_min_in & -y_min_out |
|
||||
+y_max_in & -y_max_out)
|
||||
|
||||
prism = prism & ~corners
|
||||
|
||||
return prism
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue