crp/OpenMC
1
0
Fork 0
mirror of https://github.com/openmc-dev/openmc.git synced 2026-07-28 06:05:58 -04:00
Code Issues Projects Releases Packages Wiki Activity Actions 5

Add C++ Surface pointers to Fortran Surfaces

Browse source
This commit is contained in:
Sterling Harper 2018-01-24 11:56:20 -05:00
parent 983cbd48c6
commit 75501a08e9
9 changed files with 252 additions and 203 deletions
Download patch file Download diff file Expand all files Collapse all files

68
src/geometry.F90
View file

@ -14,66 +14,6 @@ module geometry
implicit none
interface
pure function surface_sense_c(surf_ind, xyz, uvw) &
bind(C, name='surface_sense') result(sense)
use ISO_C_BINDING
implicit none
integer(C_INT), intent(in), value :: surf_ind;
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(in) :: uvw(3);
logical(C_BOOL) :: sense;
end function surface_sense_c
pure subroutine surface_reflect_c(surf_ind, xyz, uvw) &
bind(C, name='surface_reflect')
use ISO_C_BINDING
implicit none
integer(C_INT), intent(in), value :: surf_ind;
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(inout) :: uvw(3);
end subroutine surface_reflect_c
pure function surface_distance_c(surf_ind, xyz, uvw, coincident) &
bind(C, name='surface_distance') result(d)
use ISO_C_BINDING
implicit none
integer(C_INT), intent(in), value :: surf_ind;
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(in) :: uvw(3);
logical(C_BOOL), intent(in), value :: coincident;
real(C_DOUBLE) :: d;
end function surface_distance_c
pure subroutine surface_normal_c(surf_ind, xyz, uvw) &
bind(C, name='surface_normal')
use ISO_C_BINDING
implicit none
integer(C_INT), intent(in), value :: surf_ind;
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(out) :: uvw(3);
end subroutine surface_normal_c
function surface_periodic_c(surf_ind1, xyz, uvw) &
bind(C, name="surface_periodic") result(rotational)
use ISO_C_BINDING
implicit none
integer(C_INT), intent(in), value :: surf_ind1;
real(C_DOUBLE), intent(inout) :: xyz(3);
real(C_DOUBLE), intent(inout) :: uvw(3);
logical(C_BOOL) :: rotational
end function surface_periodic_c
function surface_i_periodic_c(surf_ind) bind(C, name="surface_i_periodic") &
result(i_periodic)
use ISO_C_BINDING
implicit none
integer(C_INT), intent(in), value :: surf_ind
integer(C_INT) :: i_periodic
end function
end interface
contains
!===============================================================================
@ -125,7 +65,7 @@ contains
in_cell = .false.
exit
else
actual_sense = surface_sense_c(abs(token)-1, &
actual_sense = surfaces(abs(token)) % sense( &
p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
if (actual_sense .neqv. (token > 0)) then
in_cell = .false.
@ -174,7 +114,7 @@ contains
elseif (-token == p%surface) then
stack(i_stack) = .false.
else
actual_sense = surface_sense_c(abs(token)-1, &
actual_sense = surfaces(abs(token)) % sense( &
p%coord(p%n_coord)%xyz, p%coord(p%n_coord)%uvw)
stack(i_stack) = (actual_sense .eqv. (token > 0))
end if
@ -579,7 +519,7 @@ contains
if (index_surf >= OP_UNION) cycle
! Calculate distance to surface
d = surface_distance_c(index_surf-1, p % coord(j) % xyz, &
d = surfaces(index_surf) % distance(p % coord(j) % xyz, &
p % coord(j) % uvw, logical(coincident, kind=C_BOOL))
! Check if calculated distance is new minimum
@ -799,7 +739,7 @@ contains
! traveling into if the surface is crossed
if (.not. c % simple) then
xyz_cross(:) = p % coord(j) % xyz + d_surf*p % coord(j) % uvw
call surface_normal_c(abs(level_surf_cross)-1, xyz_cross, surf_uvw)
call surfaces(abs(level_surf_cross)) % normal(xyz_cross, surf_uvw)
if (dot_product(p % coord(j) % uvw, surf_uvw) > ZERO) then
surface_crossed = abs(level_surf_cross)
else

60
src/input_xml.F90
View file

@ -924,19 +924,15 @@ contains
logical :: file_exists
logical :: boundary_exists
character(MAX_LINE_LEN) :: filename
character(MAX_WORD_LEN) :: word
character(MAX_WORD_LEN), allocatable :: sarray(:)
character(:), allocatable :: region_spec
type(Cell), pointer :: c
class(Surface), pointer :: s
class(Lattice), pointer :: lat
type(XMLDocument) :: doc
type(XMLNode) :: root
type(XMLNode) :: node_cell
type(XMLNode) :: node_surf
type(XMLNode) :: node_lat
type(XMLNode), allocatable :: node_cell_list(:)
type(XMLNode), allocatable :: node_surf_list(:)
type(XMLNode), allocatable :: node_rlat_list(:)
type(XMLNode), allocatable :: node_hlat_list(:)
type(VectorInt) :: tokens
@ -968,66 +964,18 @@ contains
! applied to a surface
boundary_exists = .false.
! get pointer to list of xml <surface>
call get_node_list(root, "surface", node_surf_list)
call read_surfaces(root % ptr)
! Get number of <surface> tags
n_surfaces = size(node_surf_list)
! Check for no surfaces
if (n_surfaces == 0) then
call fatal_error("No surfaces found in geometry.xml!")
end if
! Allocate cells array
! Allocate surfaces array
allocate(surfaces(n_surfaces))
do i = 1, n_surfaces
! Get pointer to i-th surface node
node_surf = node_surf_list(i)
surfaces(i) % ptr = surface_pointer_c(i - 1);
s => surfaces(i)
if (surfaces(i) % bc() /= BC_TRANSMIT) boundary_exists = .true.
! Copy data into cells
if (check_for_node(node_surf, "id")) then
call get_node_value(node_surf, "id", s%id)
else
call fatal_error("Must specify id of surface in geometry XML file.")
end if
! Check to make sure 'id' hasn't been used
if (surface_dict % has(s%id)) then
call fatal_error("Two or more surfaces use the same unique ID: " &
// to_str(s%id))
end if
! Check to make sure that the proper number of coefficients
! have been specified for the given type of surface. Then copy
! surface coordinates.
! Boundary conditions
word = ''
if (check_for_node(node_surf, "boundary")) &
call get_node_value(node_surf, "boundary", word)
select case (to_lower(word))
case ('transmission', 'transmit', '')
s%bc = BC_TRANSMIT
case ('vacuum')
s%bc = BC_VACUUM
boundary_exists = .true.
case ('reflective', 'reflect', 'reflecting')
s%bc = BC_REFLECT
boundary_exists = .true.
case ('periodic')
s%bc = BC_PERIODIC
boundary_exists = .true.
case default
call fatal_error("Unknown boundary condition '" // trim(word) // &
&"' specified on surface " // trim(to_str(s%id)))
end select
! Add surface to dictionary
call surface_dict % set(s%id, i)
call surface_dict % set(surfaces(i) % id(), i)
end do
! Check to make sure a boundary condition was applied to at least one

2
src/output.F90
View file

@ -259,7 +259,7 @@ contains
! Print surface
if (p % surface /= NONE) then
write(ou,*) ' Surface = ' // to_str(sign(surfaces(i)%id, p % surface))
write(ou,*) ' Surface = ' // to_str(sign(surfaces(i)%id(), p % surface))
end if
! Display weight, energy, grid index, and interpolation factor

15
src/summary.F90
View file

@ -24,17 +24,6 @@ module summary
public :: write_summary
interface
subroutine surface_to_hdf5_c(surf_ind, group) &
bind(C, name='surface_to_hdf5')
use ISO_C_BINDING
use hdf5
implicit none
integer(C_INT), intent(in), value :: surf_ind
integer(HID_T), intent(in), value :: group
end subroutine surface_to_hdf5_c
end interface
contains
!===============================================================================
@ -232,7 +221,7 @@ contains
region_spec = trim(region_spec) // " |"
case default
region_spec = trim(region_spec) // " " // to_str(&
sign(surfaces(abs(k))%id, k))
sign(surfaces(abs(k))%id(), k))
end select
end do
call write_dataset(cell_group, "region", adjustl(region_spec))
@ -250,7 +239,7 @@ contains
! Write information on each surface
SURFACE_LOOP: do i = 1, n_surfaces
call surface_to_hdf5_c(i-1, surfaces_group)
call surfaces(i) % to_hdf5(surfaces_group)
end do SURFACE_LOOP
call close_group(surfaces_group)

13
src/surface.cpp
View file

@ -1053,11 +1053,13 @@ extern "C" void
read_surfaces(pugi::xml_node *node)
{
// Count the number of surfaces.
int n_surfaces{0};
for (pugi::xml_node surf_node = node->child("surface"); surf_node;
surf_node = surf_node.next_sibling("surface")) {
n_surfaces++;
}
if (n_surfaces == 0) {
fatal_error("No surfaces found in geometry.xml!");
}
// Allocate the array of Surface pointers.
surfaces_c = new Surface* [n_surfaces];
@ -1137,7 +1139,14 @@ read_surfaces(pugi::xml_node *node)
// Downcast to the PeriodicSurface type.
Surface *surf_base = surfaces_c[i_surf];
PeriodicSurface *surf = dynamic_cast<PeriodicSurface *>(surf_base);
//TODO: check for null pointer
// Make sure this surface inherits from PeriodicSurface.
if (surf == NULL) {
std::string err_msg{"Periodic boundary condition not supported for "
"surface "};
err_msg += std::to_string(surf_base->id);
err_msg += ". Periodic BCs are only supported for planar surfaces.";
}
// See if this surface makes part of the global bounding box.
struct BoundingBox bb = surf->bounding_box();

72
src/surface.h
View file

@ -28,6 +28,8 @@ const int C_NONE{-1};
// Global variables
//==============================================================================
int n_surfaces{0};
class Surface;
Surface **surfaces_c;
@ -375,57 +377,43 @@ public:
// Fortran compatibility functions
//==============================================================================
extern "C" bool
surface_sense(int surf_ind, double xyz[3], double uvw[3])
{
return surfaces_c[surf_ind]->sense(xyz, uvw);
}
extern "C" Surface* surface_pointer(int surf_ind) {return surfaces_c[surf_ind];}
extern "C" void
surface_reflect(int surf_ind, double xyz[3], double uvw[3])
{
surfaces_c[surf_ind]->reflect(xyz, uvw);
}
extern "C" int surface_id(Surface *surf) {return surf->id;}
extern "C" int surface_bc(Surface *surf) {return surf->bc;}
extern "C" bool surface_sense(Surface *surf, double xyz[3], double uvw[3])
{return surf->sense(xyz, uvw);}
extern "C" void surface_reflect(Surface *surf, double xyz[3], double uvw[3])
{surf->reflect(xyz, uvw);}
extern "C" double
surface_distance(int surf_ind, double xyz[3], double uvw[3], bool coincident)
{
return surfaces_c[surf_ind]->distance(xyz, uvw, coincident);
}
surface_distance(Surface *surf, double xyz[3], double uvw[3], bool coincident)
{return surf->distance(xyz, uvw, coincident);}
extern "C" void
surface_normal(int surf_ind, double xyz[3], double uvw[3])
{
return surfaces_c[surf_ind]->normal(xyz, uvw);
}
extern "C" void surface_normal(Surface *surf, double xyz[3], double uvw[3])
{return surf->normal(xyz, uvw);}
extern "C" void
surface_to_hdf5(int surf_ind, hid_t group)
{
surfaces_c[surf_ind]->to_hdf5(group);
}
extern "C" void surface_to_hdf5(Surface *surf, hid_t group)
{surf->to_hdf5(group);}
extern "C" int surface_i_periodic(PeriodicSurface *surf)
{return surf->i_periodic;}
extern "C" bool
surface_periodic(int surf_ind1, double xyz[3], double uvw[3])
{
// Hopefully this function has only been called for a pair of surfaces that
// support periodic BCs (checking should have been done when reading the
// geometry XML). Downcast the surfaces to the PeriodicSurface type so we
// can call the periodic_translate method.
Surface *surf1_gen = surfaces_c[surf_ind1];
PeriodicSurface *surf1 = dynamic_cast<PeriodicSurface *>(surf1_gen);
Surface *surf2_gen = surfaces_c[surf1->i_periodic];
PeriodicSurface *surf2 = dynamic_cast<PeriodicSurface *>(surf2_gen);
surface_periodic(PeriodicSurface *surf, PeriodicSurface *other, double xyz[3],
double uvw[3])
{return surf->periodic_translate(other, xyz, uvw);}
// Call the type-bound methods.
return surf2->periodic_translate(surf1, xyz, uvw);
}
extern "C" int
surface_i_periodic(int surf_ind)
extern "C" void free_memory_surfaces_c()
{
PeriodicSurface *surf = dynamic_cast<PeriodicSurface *>(surfaces_c[surf_ind]);
return surf->i_periodic;
for (int i = 0; i < n_surfaces; i++) {delete surfaces_c[i];}
delete surfaces_c;
surfaces_c = NULL;
n_surfaces = 0;
surface_dict.clear();
}
#endif // SURFACE_H

180
src/surface_header.F90
View file

@ -1,26 +1,132 @@
module surface_header
use, intrinsic :: ISO_C_BINDING
use hdf5
use constants, only: NONE
use dict_header, only: DictIntInt
implicit none
interface
pure function surface_pointer_c(surf_ind) &
bind(C, name='surface_pointer') result(ptr)
use ISO_C_BINDING
implicit none
integer(C_INT), intent(in), value :: surf_ind
type(C_PTR) :: ptr
end function surface_pointer_c
pure function surface_id_c(surf_ptr) bind(C, name='surface_id') result(id)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
integer(C_INT) :: id
end function surface_id_c
pure function surface_bc_c(surf_ptr) bind(C, name='surface_bc') result(bc)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
integer(C_INT) :: bc
end function surface_bc_c
pure function surface_sense_c(surf_ptr, xyz, uvw) &
bind(C, name='surface_sense') result(sense)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
real(C_DOUBLE), intent(in) :: xyz(3)
real(C_DOUBLE), intent(in) :: uvw(3)
logical(C_BOOL) :: sense
end function surface_sense_c
pure subroutine surface_reflect_c(surf_ptr, xyz, uvw) &
bind(C, name='surface_reflect')
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(inout) :: uvw(3);
end subroutine surface_reflect_c
pure function surface_distance_c(surf_ptr, xyz, uvw, coincident) &
bind(C, name='surface_distance') result(d)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(in) :: uvw(3);
logical(C_BOOL), intent(in), value :: coincident;
real(C_DOUBLE) :: d;
end function surface_distance_c
pure subroutine surface_normal_c(surf_ptr, xyz, uvw) &
bind(C, name='surface_normal')
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(out) :: uvw(3);
end subroutine surface_normal_c
subroutine surface_to_hdf5_c(surf_ptr, group) &
bind(C, name='surface_to_hdf5')
use ISO_C_BINDING
use hdf5
implicit none
type(C_PTR), intent(in), value :: surf_ptr
integer(HID_T), intent(in), value :: group
end subroutine surface_to_hdf5_c
pure function surface_i_periodic_c(surf_ptr) &
bind(C, name="surface_i_periodic") result(i_periodic)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
integer(C_INT) :: i_periodic
end function surface_i_periodic_c
function surface_periodic_c(surf_ptr, other_ptr, xyz, uvw) &
bind(C, name="surface_periodic") result(rotational)
use ISO_C_BINDING
implicit none
type(C_PTR), intent(in), value :: surf_ptr
type(C_PTR), intent(in), value :: other_ptr
real(C_DOUBLE), intent(inout) :: xyz(3);
real(C_DOUBLE), intent(inout) :: uvw(3);
logical(C_BOOL) :: rotational
end function surface_periodic_c
subroutine free_memory_surfaces_c() bind(C)
end subroutine free_memory_surfaces_c
end interface
!===============================================================================
! SURFACE type defines a first- or second-order surface that can be used to
! construct closed volumes (cells)
!===============================================================================
type :: Surface
integer :: id ! Unique ID
integer, allocatable :: &
neighbor_pos(:), & ! List of cells on positive side
neighbor_neg(:) ! List of cells on negative side
integer :: bc ! Boundary condition
type(C_PTR) :: ptr
contains
procedure :: id => surface_id
procedure :: bc => surface_bc
procedure :: sense => surface_sense
procedure :: reflect => surface_reflect
procedure :: distance => surface_distance
procedure :: normal => surface_normal
procedure :: to_hdf5 => surface_to_hdf5
procedure :: i_periodic => surface_i_periodic
procedure :: periodic_translate => surface_periodic
end type Surface
integer(C_INT32_T), bind(C) :: n_surfaces ! # of surfaces
integer(C_INT), bind(C) :: n_surfaces ! # of surfaces
type(Surface), allocatable, target :: surfaces(:)
@ -29,14 +135,78 @@ module surface_header
contains
pure function surface_id(this) result(id)
class(Surface), intent(in) :: this
integer(C_INT) :: id
id = surface_id_c(this % ptr)
end function surface_id
pure function surface_bc(this) result(bc)
class(Surface), intent(in) :: this
integer(C_INT) :: bc
bc = surface_bc_c(this % ptr)
end function surface_bc
pure function surface_sense(this, xyz, uvw) result(sense)
class(Surface), intent(in) :: this
real(C_DOUBLE), intent(in) :: xyz(3)
real(C_DOUBLE), intent(in) :: uvw(3)
logical(C_BOOL) :: sense
sense = surface_sense_c(this % ptr, xyz, uvw)
end function surface_sense
pure subroutine surface_reflect(this, xyz, uvw)
class(Surface), intent(in) :: this
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(inout) :: uvw(3);
call surface_reflect_c(this % ptr, xyz, uvw)
end subroutine surface_reflect
pure function surface_distance(this, xyz, uvw, coincident) result(d)
class(Surface), intent(in) :: this
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(in) :: uvw(3);
logical(C_BOOL), intent(in) :: coincident;
real(C_DOUBLE) :: d;
d = surface_distance_c(this % ptr, xyz, uvw, coincident)
end function surface_distance
pure subroutine surface_normal(this, xyz, uvw)
class(Surface), intent(in) :: this
real(C_DOUBLE), intent(in) :: xyz(3);
real(C_DOUBLE), intent(out) :: uvw(3);
call surface_normal_c(this % ptr, xyz, uvw)
end subroutine surface_normal
subroutine surface_to_hdf5(this, group)
class(Surface), intent(in) :: this
integer(HID_T), intent(in) :: group
call surface_to_hdf5_c(this % ptr, group)
end subroutine surface_to_hdf5
pure function surface_i_periodic(this) result(i_periodic)
class(Surface), intent(in) :: this
integer(C_INT) :: i_periodic
i_periodic = surface_i_periodic_c(this % ptr)
end function surface_i_periodic
function surface_periodic(this, other, xyz, uvw) result(rotational)
class(Surface), intent(in) :: this
class(Surface), intent(in) :: other
real(C_DOUBLE), intent(inout) :: xyz(3);
real(C_DOUBLE), intent(inout) :: uvw(3);
logical(C_BOOL) :: rotational
rotational = surface_periodic_c(this % ptr, other % ptr, xyz, uvw)
end function surface_periodic
!===============================================================================
! FREE_MEMORY_SURFACES deallocates global arrays defined in this module
!===============================================================================
subroutine free_memory_surfaces()
n_surfaces = 0
if (allocated(surfaces)) deallocate(surfaces)
call surface_dict % clear()
call free_memory_surfaces_c()
end subroutine free_memory_surfaces
end module surface_header

2
src/tallies/tally_filter_surface.F90
View file

@ -105,7 +105,7 @@ contains
integer, intent(in) :: bin
character(MAX_LINE_LEN) :: label
label = "Surface " // to_str(surfaces(this % surfaces(bin)) % id)
label = "Surface " // to_str(surfaces(this % surfaces(bin)) % id())
end function text_label_surface
end module tally_filter_surface

43
src/tracking.F90
View file

@ -4,9 +4,8 @@ module tracking
use cross_section, only: calculate_xs
use error, only: fatal_error, warning, write_message
use geometry_header, only: cells
use geometry, only: find_cell, distance_to_boundary, cross_lattice, &
check_cell_overlap, surface_reflect_c, &
surface_periodic_c, surface_i_periodic_c
use geometry, only: find_cell, distance_to_boundary, cross_lattice,&
check_cell_overlap
use message_passing
use mgxs_header
use nuclide_header
@ -296,14 +295,15 @@ contains
logical :: rotational ! if rotational periodic BC applied
logical :: found ! particle found in universe?
class(Surface), pointer :: surf
class(Surface), pointer :: surf2 ! periodic partner surface
i_surface = abs(p % surface)
surf => surfaces(i_surface)
if (verbosity >= 10 .or. trace) then
call write_message(" Crossing surface " // trim(to_str(surf % id)))
call write_message(" Crossing surface " // trim(to_str(surf % id())))
end if
if (surf % bc == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then
if (surf % bc() == BC_VACUUM .and. (run_mode /= MODE_PLOTTING)) then
! =======================================================================
! PARTICLE LEAKS OUT OF PROBLEM
@ -328,11 +328,11 @@ contains
! Display message
if (verbosity >= 10 .or. trace) then
call write_message(" Leaked out of surface " &
&// trim(to_str(surf % id)))
&// trim(to_str(surf % id())))
end if
return
elseif (surf % bc == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then
elseif (surf % bc() == BC_REFLECT .and. (run_mode /= MODE_PLOTTING)) then
! =======================================================================
! PARTICLE REFLECTS FROM SURFACE
@ -355,7 +355,7 @@ contains
end if
! Reflect particle off surface
call surface_reflect_c(i_surface-1, p%coord(1)%xyz, p%coord(1)%uvw)
call surf % reflect(p%coord(1)%xyz, p%coord(1)%uvw)
! Make sure new particle direction is normalized
u = p%coord(1)%uvw(1)
@ -376,7 +376,7 @@ contains
call find_cell(p, found)
if (.not. found) then
call p % mark_as_lost("Couldn't find particle after reflecting&
& from surface " // trim(to_str(surf % id)) // ".")
& from surface " // trim(to_str(surf % id())) // ".")
return
end if
@ -386,10 +386,10 @@ contains
! Diagnostic message
if (verbosity >= 10 .or. trace) then
call write_message(" Reflected from surface " &
&// trim(to_str(surf%id)))
&// trim(to_str(surf%id())))
end if
return
elseif (surf % bc == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then
elseif (surf % bc() == BC_PERIODIC .and. run_mode /= MODE_PLOTTING) then
! =======================================================================
! PERIODIC BOUNDARY
@ -404,7 +404,6 @@ contains
! Score surface currents since reflection causes the direction of the
! particle to change -- artificially move the particle slightly back in
! case the surface crossing is coincident with a mesh boundary
if (active_current_tallies % size() > 0) then
xyz = p % coord(1) % xyz
p % coord(1) % xyz = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw
@ -412,14 +411,19 @@ contains
p % coord(1) % xyz = xyz
end if
rotational = surface_periodic_c(i_surface-1, p % coord(1) % xyz, &
p % coord(1) % uvw)
! Get a pointer to the partner periodic surface. Offset the index to
! correct for C vs. Fortran indexing.
surf2 => surfaces(surf % i_periodic() + 1)
! Adjust the particle's location and direction.
rotational = surf2 % periodic_translate(surf, p % coord(1) % xyz, &
p % coord(1) % uvw)
! Reassign particle's surface
if (rotational) then
p % surface = surface_i_periodic_c(i_surface-1) + 1
p % surface = surf % i_periodic() + 1
else
p % surface = sign(surface_i_periodic_c(i_surface-1) + 1, p % surface)
p % surface = sign(surf % i_periodic() + 1, p % surface)
end if
! Figure out what cell particle is in now
@ -427,7 +431,8 @@ contains
call find_cell(p, found)
if (.not. found) then
call p % mark_as_lost("Couldn't find particle after hitting &
&periodic boundary on surface " // trim(to_str(surf % id)) // ".")
&periodic boundary on surface " // trim(to_str(surf % id())) &
// ".")
return
end if
@ -437,7 +442,7 @@ contains
! Diagnostic message
if (verbosity >= 10 .or. trace) then
call write_message(" Hit periodic boundary on surface " &
// trim(to_str(surf%id)))
// trim(to_str(surf%id())))
end if
return
end if
@ -484,7 +489,7 @@ contains
if (.not. found) then
call p % mark_as_lost("After particle " // trim(to_str(p % id)) &
// " crossed surface " // trim(to_str(surf % id)) &
// " crossed surface " // trim(to_str(surf % id())) &
// " it could not be located in any cell and it did not leak.")
return
end if