mirror of
https://github.com/openmc-dev/openmc.git
synced 2026-07-28 14:15:42 -04:00
Merge pull request #879 from paulromano/prism-fixes
Bugfix for applying bounding type in get_hexagonal_prism
This commit is contained in:
sam
GitHub
commit
d17942de68
3 changed files with 136 additions and 265 deletions
|
|
@ -1,6 +1,8 @@
|
|||
from __future__ import division
|
||||
from abc import ABCMeta
|
||||
from collections import Iterable, 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
|
||||
|
|
@ -1326,6 +1328,7 @@ class Sphere(Surface):
|
|||
z = point[2] - self.z0
|
||||
return x**2 + y**2 + z**2 - self.r**2
|
||||
|
||||
|
||||
@add_metaclass(ABCMeta)
|
||||
class Cone(Surface):
|
||||
"""A conical surface parallel to the x-, y-, or z-axis.
|
||||
|
|
@ -1995,147 +1998,66 @@ def get_rectangular_prism(width, height, axis='z', origin=(0., 0.),
|
|||
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_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':
|
||||
min_y = YPlane(name='minimum y', y0=-width/2.+origin[0],
|
||||
boundary_type=boundary_type)
|
||||
max_y = YPlane(name='maximum y', y0=+width/2.+origin[0],
|
||||
boundary_type=boundary_type)
|
||||
min_z = ZPlane(name='minimum z', z0=-height/2.+origin[1],
|
||||
boundary_type=boundary_type)
|
||||
max_z = ZPlane(name='maximum z', z0=+height/2.+origin[1],
|
||||
boundary_type=boundary_type)
|
||||
prism = +min_y & -max_y & +min_z & -max_z
|
||||
x1, x2 = 'y', 'z'
|
||||
elif axis == 'y':
|
||||
min_x = XPlane(name='minimum x', x0=-width/2.+origin[0],
|
||||
boundary_type=boundary_type)
|
||||
max_x = XPlane(name='maximum x', x0=+width/2.+origin[0],
|
||||
boundary_type=boundary_type)
|
||||
min_z = ZPlane(name='minimum z', z0=-height/2.+origin[1],
|
||||
boundary_type=boundary_type)
|
||||
max_z = ZPlane(name='maximum z', z0=+height/2.+origin[1],
|
||||
boundary_type=boundary_type)
|
||||
prism = +min_x & -max_x & +min_z & -max_z
|
||||
x1, x2 = 'x', 'z'
|
||||
else:
|
||||
min_x = XPlane(name='minimum x', x0=-width/2.+origin[0],
|
||||
boundary_type=boundary_type)
|
||||
max_x = XPlane(name='maximum x', x0=+width/2.+origin[0],
|
||||
boundary_type=boundary_type)
|
||||
min_y = YPlane(name='minimum y', y0=-height/2.+origin[1],
|
||||
boundary_type=boundary_type)
|
||||
max_y = YPlane(name='maximum y', y0=+height/2.+origin[1],
|
||||
boundary_type=boundary_type)
|
||||
prism = +min_x & -max_x & +min_y & -max_y
|
||||
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.:
|
||||
if axis == 'x':
|
||||
y_min_z_min = XCylinder(name='y min z min',
|
||||
y0=origin[0] - width /2. + corner_radius,
|
||||
z0=origin[1] - height/2. + corner_radius,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
y_min_z_max = XCylinder(name='y min z max',
|
||||
y0=origin[0] - width /2. + corner_radius,
|
||||
z0=origin[1] + height/2. - corner_radius,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
y_max_z_min = XCylinder(name='y max z min',
|
||||
y0=origin[0] + width /2. - corner_radius,
|
||||
z0=origin[1] - height/2. + corner_radius,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
y_max_z_max = XCylinder(name='y max z max',
|
||||
y0=origin[0] + width /2. - corner_radius,
|
||||
z0=origin[1] + height/2. - corner_radius,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
y_min = YPlane(name='min y', y0=-width/2.+origin[0]+corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
y_max = YPlane(name='max y', y0=+width/2.+origin[0]-corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
z_min = ZPlane(name='min z', z0=-height/2.+origin[1]+corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
z_max = ZPlane(name='max z', z0=+height/2.+origin[1]-corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
args = {'R': corner_radius, 'boundary_type': boundary_type}
|
||||
|
||||
corners = (+y_min_z_min & -y_min & -z_min) | \
|
||||
(+y_min_z_max & -y_min & +z_max) | \
|
||||
(+y_max_z_min & +y_max & -z_min) | \
|
||||
(+y_max_z_max & +y_max & +z_max)
|
||||
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)
|
||||
|
||||
elif axis == 'y':
|
||||
x_min_z_min = YCylinder(name='x min z min',
|
||||
x0=origin[0] - width /2. + corner_radius,
|
||||
z0=origin[1] - height/2. + corner_radius,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_min_z_max = YCylinder(name='x min z max',
|
||||
x0=origin[0] - width /2. + corner_radius,
|
||||
z0=origin[1] + height/2. - corner_radius,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max_z_min = YCylinder(name='x max z min',
|
||||
x0=origin[0] + width /2. - corner_radius,
|
||||
z0=origin[1] - height/2. + corner_radius,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max_z_max = YCylinder(name='x max z max',
|
||||
x0=origin[0] + width /2. - corner_radius,
|
||||
z0=origin[1] + height/2. - corner_radius,
|
||||
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)
|
||||
|
||||
x_min = XPlane(name='min x', x0=-width/2.+origin[0]+corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max = XPlane(name='max x', x0=+width/2.+origin[0]-corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
z_min = ZPlane(name='min z', z0=-height/2.+origin[1]+corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
z_max = ZPlane(name='max z', z0=+height/2.+origin[1]-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_max = cyl(name='{} min {} max'.format(x1, x2), **args)
|
||||
|
||||
corners = (+x_min_z_min & -x_min & -z_min) | \
|
||||
(+x_min_z_max & -x_min & +z_max) | \
|
||||
(+x_max_z_min & +x_max & -z_min) | \
|
||||
(+x_max_z_max & +x_max & +z_max)
|
||||
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)
|
||||
|
||||
else:
|
||||
x_min_y_min = ZCylinder(name='x min y min',
|
||||
x0=origin[0] - width /2. + corner_radius,
|
||||
y0=origin[1] - height/2. + corner_radius,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_min_y_max = ZCylinder(name='x min y max',
|
||||
x0=origin[0] - width /2. + corner_radius,
|
||||
y0=origin[1] + height/2. - corner_radius,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max_y_min = ZCylinder(name='x max y min',
|
||||
x0=origin[0] + width /2. - corner_radius,
|
||||
y0=origin[1] - height/2. + corner_radius,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max_y_max = ZCylinder(name='x max y max',
|
||||
x0=origin[0] + width /2. - corner_radius,
|
||||
y0=origin[1] + height/2. - corner_radius,
|
||||
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_max_x2_max = cyl(name='{} max {} max'.format(x1, x2), **args)
|
||||
|
||||
x_min = XPlane(name='min x', x0=-width/2.+origin[0]+corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max = XPlane(name='max x', x0=+width/2.+origin[0]-corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
y_min = YPlane(name='min y', y0=-height/2.+origin[1]+corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
y_max = YPlane(name='max y', y0=+height/2.+origin[1]-corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
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 = (+x_min_y_min & -x_min & -y_min) | \
|
||||
(+x_min_y_max & -x_min & +y_max) | \
|
||||
(+x_max_y_min & +x_max & -y_min) | \
|
||||
(+x_max_y_max & +x_max & +y_max)
|
||||
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
|
||||
|
||||
|
|
@ -2170,8 +2092,8 @@ def get_hexagonal_prism(edge_length=1., orientation='y',
|
|||
l = edge_length
|
||||
|
||||
if orientation == 'y':
|
||||
right = XPlane(x0=sqrt(3.)/2.*l)
|
||||
left = XPlane(x0=-sqrt(3.)/2.*l)
|
||||
right = XPlane(x0=sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
left = XPlane(x0=-sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
c = sqrt(3.)/3.
|
||||
|
||||
# y = -x/sqrt(3) + a
|
||||
|
|
@ -2188,9 +2110,14 @@ def get_hexagonal_prism(edge_length=1., orientation='y',
|
|||
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=sqrt(3.)/2.*l)
|
||||
bottom = YPlane(y0=-sqrt(3.)/2.*l)
|
||||
top = YPlane(y0=sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
bottom = YPlane(y0=-sqrt(3.)/2*l, boundary_type=boundary_type)
|
||||
c = sqrt(3.)
|
||||
|
||||
# y = -sqrt(3)*(x - a)
|
||||
|
|
@ -2207,148 +2134,68 @@ def get_hexagonal_prism(edge_length=1., orientation='y',
|
|||
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':
|
||||
c = sqrt(3.)
|
||||
x_min_y_min_in = ZCylinder(name='x min y min in',
|
||||
x0=-0.5 * (l - 2./c*corner_radius),
|
||||
y0=-c/2. * (l - 2./c*corner_radius),
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_min_y_max_in = ZCylinder(name='x min y max in',
|
||||
x0= 0.5 * (l - 2./c*corner_radius),
|
||||
y0=-c/2. * (l - 2./c*corner_radius),
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max_y_min_in = ZCylinder(name='x max y min in',
|
||||
x0=-0.5 * (l - 2./c*corner_radius),
|
||||
y0= c/2. * (l - 2./c*corner_radius),
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max_y_max_in = ZCylinder(name='x max y max in',
|
||||
x0= 0.5 * (l - 2./c*corner_radius),
|
||||
y0= c/2. * (l - 2./c*corner_radius),
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_min_in = ZCylinder(name='x min in',
|
||||
x0=-(l - 2./c*corner_radius),
|
||||
y0=0.,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max_in = ZCylinder(name='x max in',
|
||||
x0= (l - 2./c*corner_radius),
|
||||
y0=0.,
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_min_y_min_in = cyl1(name='x min y min in', x0=-t/2, y0=-c*t)
|
||||
x_min_y_max_in = cyl1(name='x min y max in', x0=t/2, y0=-c*t)
|
||||
x_max_y_min_in = cyl1(name='x max y min in', x0=-t/2, y0=c*t)
|
||||
x_max_y_max_in = cyl1(name='x max y max in', x0=t/2, y0=c*t)
|
||||
x_min_in = cyl1(name='x min in', x0=-t)
|
||||
x_max_in = cyl1(name='x max in', x0=t)
|
||||
|
||||
x_min_y_min_out = ZCylinder(name='x min y min out',
|
||||
x0=-0.5 * l,
|
||||
y0=-c/2. * l,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
x_min_y_max_out = ZCylinder(name='x min y max out',
|
||||
x0=+0.5 * l,
|
||||
y0=-c/2. * l,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
x_max_y_min_out = ZCylinder(name='x max y min out',
|
||||
x0=-0.5 * l,
|
||||
y0= c/2. * l,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
x_max_y_max_out = ZCylinder(name='x max y max out',
|
||||
x0= 0.5 * l,
|
||||
y0= c/2. * l,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
x_min_out = ZCylinder(name='x min out',
|
||||
x0=-l,
|
||||
y0=0.,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
x_max_out = ZCylinder(name='x max out',
|
||||
x0= l,
|
||||
y0=0.,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
x_min_y_min_out = cyl2(name='x min y min out', x0=-l/2, y0=-c*l)
|
||||
x_min_y_max_out = cyl2(name='x min y max out', x0=l/2, y0=-c*l)
|
||||
x_max_y_min_out = cyl2(name='x max y min out', x0=-l/2, y0=c*l)
|
||||
x_max_y_max_out = cyl2(name='x max y max out', x0=l/2, y0=c*l)
|
||||
x_min_out = cyl2(name='x min out', x0=-l)
|
||||
x_max_out = cyl2(name='x max out', x0=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) | \
|
||||
(+x_min_in & -x_min_out ) | \
|
||||
(+x_max_in & -x_max_out )
|
||||
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':
|
||||
c = sqrt(3.)
|
||||
x_min_y_min_in = ZCylinder(name='x min y min in',
|
||||
x0=-c/2. * (l - 2./c*corner_radius),
|
||||
y0=-0.5 * (l - 2./c*corner_radius),
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_min_y_max_in = ZCylinder(name='x min y max in',
|
||||
x0=-c/2. * (l - 2./c*corner_radius),
|
||||
y0= 0.5 * (l - 2./c*corner_radius),
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max_y_min_in = ZCylinder(name='x max y min in',
|
||||
x0= c/2. * (l - 2./c*corner_radius),
|
||||
y0=-0.5 * (l - 2./c*corner_radius),
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_max_y_max_in = ZCylinder(name='x max y max in',
|
||||
x0= c/2. * (l - 2./c*corner_radius),
|
||||
y0= 0.5 * (l - 2./c*corner_radius),
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
y_min_in = ZCylinder(name='y min in',
|
||||
x0=0.,
|
||||
y0=-(l - 2./c*corner_radius),
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
y_max_in = ZCylinder(name='y max in',
|
||||
x0=0.,
|
||||
y0= (l - 2./c*corner_radius),
|
||||
R=corner_radius,
|
||||
boundary_type=boundary_type)
|
||||
x_min_y_min_in = cyl1(name='x min y min in', x0=-c*t, y0=-t/2)
|
||||
x_min_y_max_in = cyl1(name='x min y max in', x0=-c*t, y0=t/2)
|
||||
x_max_y_min_in = cyl1(name='x max y min in', x0=c*t, y0=-t/2)
|
||||
x_max_y_max_in = cyl1(name='x max y max in', x0=c*t, y0=t/2)
|
||||
y_min_in = cyl1(name='y min in', y0=-t)
|
||||
y_max_in = cyl1(name='y max in', y0=t)
|
||||
|
||||
x_min_y_min_out = ZCylinder(name='x min y min out',
|
||||
x0=-c/2. * l,
|
||||
y0=-0.5 * l,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
x_min_y_max_out = ZCylinder(name='x min y max out',
|
||||
x0=-c/2. * l,
|
||||
y0=+0.5 * l,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
x_max_y_min_out = ZCylinder(name='x max y min out',
|
||||
x0= c/2. * l,
|
||||
y0=-0.5 * l,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
x_max_y_max_out = ZCylinder(name='x max y max out',
|
||||
x0= c/2. * l,
|
||||
y0= 0.5 * l,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
y_min_out = ZCylinder(name='y min out',
|
||||
x0=0.,
|
||||
y0=-l,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
y_max_out = ZCylinder(name='y max out',
|
||||
x0=0.,
|
||||
y0= l,
|
||||
R=corner_radius/c,
|
||||
boundary_type=boundary_type)
|
||||
x_min_y_min_out = cyl2(name='x min y min out', x0=-c*l, y0=-l/2)
|
||||
x_min_y_max_out = cyl2(name='x min y max out', x0=-c*l, y0=l/2)
|
||||
x_max_y_min_out = cyl2(name='x max y min out', x0=c*l, y0=-l/2)
|
||||
x_max_y_max_out = cyl2(name='x max y max out', x0=c*l, y0=l/2)
|
||||
y_min_out = cyl2(name='y min out', y0=-l)
|
||||
y_max_out = cyl2(name='y max out', y0=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 )
|
||||
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
|
||||
|
||||
|
|
|
|||
|
|
@ -392,6 +392,7 @@ contains
|
|||
real(8) :: v ! y-component of direction
|
||||
real(8) :: w ! z-component of direction
|
||||
real(8) :: norm ! "norm" of surface normal
|
||||
real(8) :: d ! distance between point and plane
|
||||
real(8) :: xyz(3) ! Saved global coordinate
|
||||
integer :: i_surface ! index in surfaces
|
||||
logical :: found ! particle found in universe?
|
||||
|
|
@ -532,6 +533,20 @@ contains
|
|||
type is (SurfaceZPlane)
|
||||
p % coord(1) % xyz(3) = opposite % z0
|
||||
end select
|
||||
|
||||
type is (SurfacePlane)
|
||||
select type (opposite => surfaces(surf % i_periodic) % obj)
|
||||
type is (SurfacePlane)
|
||||
! Get surface normal for opposite plane
|
||||
xyz(:) = opposite % normal(p % coord(1) % xyz)
|
||||
|
||||
! Determine distance to plane
|
||||
norm = xyz(1)*xyz(1) + xyz(2)*xyz(2) + xyz(3)*xyz(3)
|
||||
d = opposite % evaluate(p % coord(1) % xyz) / norm
|
||||
|
||||
! Move particle along normal vector based on distance
|
||||
p % coord(1) % xyz(:) = p % coord(1) % xyz(:) - d*xyz
|
||||
end select
|
||||
end select
|
||||
|
||||
! Reassign particle's surface
|
||||
|
|
|
|||
|
|
@ -1658,6 +1658,15 @@ contains
|
|||
surf % i_periodic = surface_dict % get_key(surf % i_periodic)
|
||||
end if
|
||||
|
||||
type is (SurfacePlane)
|
||||
if (surf % i_periodic == NONE) then
|
||||
call fatal_error("No matching periodic surface specified for &
|
||||
&periodic boundary condition on surface " // &
|
||||
trim(to_str(surf % id)) // ".")
|
||||
else
|
||||
surf % i_periodic = surface_dict % get_key(surf % i_periodic)
|
||||
end if
|
||||
|
||||
class default
|
||||
call fatal_error("Periodic boundary condition applied to &
|
||||
&non-planar surface.")
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue