diff --git a/openmc/surface.py b/openmc/surface.py index cb5efe4302..5a8dc53564 100644 --- a/openmc/surface.py +++ b/openmc/surface.py @@ -2110,6 +2110,11 @@ 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, boundary_type=boundary_type) bottom = YPlane(y0=-sqrt(3.)/2*l, boundary_type=boundary_type) @@ -2129,6 +2134,11 @@ 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.: c = sqrt(3.)/2 diff --git a/src/geometry.F90 b/src/geometry.F90 index 4856159913..5276edb24c 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 6c1a40ee3a..b05d560752 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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.")