diff --git a/src/geometry.F90 b/src/geometry.F90 index d1519f8d21..0bacd7b5d1 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -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 diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 9b6e1c6fce..55eeca9e14 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -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 - call get_node_list(root, "surface", node_surf_list) call read_surfaces(root % ptr) - ! Get number of 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 diff --git a/src/output.F90 b/src/output.F90 index 4113eea472..316fe95941 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -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 diff --git a/src/summary.F90 b/src/summary.F90 index 9b82a03967..3aeb42178b 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -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) diff --git a/src/surface.cpp b/src/surface.cpp index ed111b0370..369334d500 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -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(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(); diff --git a/src/surface.h b/src/surface.h index 22eb798dcb..469d884862 100644 --- a/src/surface.h +++ b/src/surface.h @@ -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(surf1_gen); - Surface *surf2_gen = surfaces_c[surf1->i_periodic]; - PeriodicSurface *surf2 = dynamic_cast(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(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 diff --git a/src/surface_header.F90 b/src/surface_header.F90 index ffe0a6db41..8985e8406e 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -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 diff --git a/src/tallies/tally_filter_surface.F90 b/src/tallies/tally_filter_surface.F90 index 6a2afa3ec6..df3bf36031 100644 --- a/src/tallies/tally_filter_surface.F90 +++ b/src/tallies/tally_filter_surface.F90 @@ -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 diff --git a/src/tracking.F90 b/src/tracking.F90 index 98bca2a279..75555c7974 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -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