diff --git a/src/geometry.F90 b/src/geometry.F90 index f1feb3ba30..d1519f8d21 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -54,16 +54,24 @@ module geometry real(C_DOUBLE), intent(out) :: uvw(3); end subroutine surface_normal_c - function surface_periodic_c(surf_ind1, surf_ind2, xyz, uvw) & + 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; - integer(C_INT), intent(in), value :: surf_ind2; 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 @@ -857,15 +865,15 @@ contains ! Copy positive neighbors to Surface instance n = neighbor_pos(i)%size() if (n > 0) then - allocate(surfaces(i)%obj%neighbor_pos(n)) - surfaces(i)%obj%neighbor_pos(:) = neighbor_pos(i)%data(1:n) + allocate(surfaces(i)%neighbor_pos(n)) + surfaces(i)%neighbor_pos(:) = neighbor_pos(i)%data(1:n) end if ! Copy negative neighbors to Surface instance n = neighbor_neg(i)%size() if (n > 0) then - allocate(surfaces(i)%obj%neighbor_neg(n)) - surfaces(i)%obj%neighbor_neg(:) = neighbor_neg(i)%data(1:n) + allocate(surfaces(i)%neighbor_neg(n)) + surfaces(i)%neighbor_neg(:) = neighbor_neg(i)%data(1:n) end if end do diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 3f22522c2a..9b6e1c6fce 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -919,12 +919,8 @@ contains integer :: id integer :: univ_id integer :: n_cells_in_univ - integer :: coeffs_reqd - integer :: i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax - real(8) :: xmin, xmax, ymin, ymax, zmin, zmax integer, allocatable :: temp_int_array(:) real(8) :: phi, theta, psi - real(8), allocatable :: coeffs(:) logical :: file_exists logical :: boundary_exists character(MAX_LINE_LEN) :: filename @@ -984,13 +980,6 @@ contains call fatal_error("No surfaces found in geometry.xml!") end if - xmin = INFINITY - xmax = -INFINITY - ymin = INFINITY - ymax = -INFINITY - zmin = INFINITY - zmax = -INFINITY - ! Allocate cells array allocate(surfaces(n_surfaces)) @@ -998,52 +987,7 @@ contains ! Get pointer to i-th surface node node_surf = node_surf_list(i) - ! Copy and interpret surface type - word = '' - if (check_for_node(node_surf, "type")) & - call get_node_value(node_surf, "type", word) - select case(to_lower(word)) - case ('x-plane') - coeffs_reqd = 1 - allocate(SurfaceXPlane :: surfaces(i)%obj) - case ('y-plane') - coeffs_reqd = 1 - allocate(SurfaceYPlane :: surfaces(i)%obj) - case ('z-plane') - coeffs_reqd = 1 - allocate(SurfaceZPlane :: surfaces(i)%obj) - case ('plane') - coeffs_reqd = 4 - allocate(SurfacePlane :: surfaces(i)%obj) - case ('x-cylinder') - coeffs_reqd = 3 - allocate(SurfaceXCylinder :: surfaces(i)%obj) - case ('y-cylinder') - coeffs_reqd = 3 - allocate(SurfaceYCylinder :: surfaces(i)%obj) - case ('z-cylinder') - coeffs_reqd = 3 - allocate(SurfaceZCylinder :: surfaces(i)%obj) - case ('sphere') - coeffs_reqd = 4 - allocate(SurfaceSphere :: surfaces(i)%obj) - case ('x-cone') - coeffs_reqd = 4 - allocate(SurfaceXCone :: surfaces(i)%obj) - case ('y-cone') - coeffs_reqd = 4 - allocate(SurfaceYCone :: surfaces(i)%obj) - case ('z-cone') - coeffs_reqd = 4 - allocate(SurfaceZCone :: surfaces(i)%obj) - case ('quadric') - coeffs_reqd = 10 - allocate(SurfaceQuadric :: surfaces(i)%obj) - case default - call fatal_error("Invalid surface type: " // trim(word)) - end select - - s => surfaces(i)%obj + s => surfaces(i) ! Copy data into cells if (check_for_node(node_surf, "id")) then @@ -1058,44 +1002,10 @@ contains // to_str(s%id)) end if - ! Copy surface name - if (check_for_node(node_surf, "name")) then - call get_node_value(node_surf, "name", s%name) - 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. - n = node_word_count(node_surf, "coeffs") - - allocate(coeffs(n)) - call get_node_array(node_surf, "coeffs", coeffs) - - select type(s) - type is (SurfaceXPlane) - ! Determine outer surfaces - xmin = min(xmin, coeffs(1)) - xmax = max(xmax, coeffs(1)) - if (xmin == coeffs(1)) i_xmin = i - if (xmax == coeffs(1)) i_xmax = i - type is (SurfaceYPlane) - ! Determine outer surfaces - ymin = min(ymin, coeffs(1)) - ymax = max(ymax, coeffs(1)) - if (ymin == coeffs(1)) i_ymin = i - if (ymax == coeffs(1)) i_ymax = i - type is (SurfaceZPlane) - ! Determine outer surfaces - zmin = min(zmin, coeffs(1)) - zmax = max(zmax, coeffs(1)) - if (zmin == coeffs(1)) i_zmin = i - if (zmax == coeffs(1)) i_zmax = i - end select - - ! No longer need coefficients - deallocate(coeffs) - ! Boundary conditions word = '' if (check_for_node(node_surf, "boundary")) & @@ -1112,12 +1022,6 @@ contains case ('periodic') s%bc = BC_PERIODIC boundary_exists = .true. - - ! Check for specification of periodic surface - if (check_for_node(node_surf, "periodic_surface_id")) then - call get_node_value(node_surf, "periodic_surface_id", & - s % i_periodic) - end if case default call fatal_error("Unknown boundary condition '" // trim(word) // & &"' specified on surface " // trim(to_str(s%id))) @@ -1134,76 +1038,6 @@ contains end if end if - ! Determine opposite side for periodic boundaries - do i = 1, size(surfaces) - if (surfaces(i) % obj % bc == BC_PERIODIC) then - select type (surf => surfaces(i) % obj) - type is (SurfaceXPlane) - if (surf % i_periodic == NONE) then - if (i == i_xmin) then - surf % i_periodic = i_xmax - elseif (i == i_xmax) then - surf % i_periodic = i_xmin - else - call fatal_error("Periodic boundary condition applied to & - &interior surface.") - end if - else - surf % i_periodic = surface_dict % get(surf % i_periodic) - end if - - type is (SurfaceYPlane) - if (surf % i_periodic == NONE) then - if (i == i_ymin) then - surf % i_periodic = i_ymax - elseif (i == i_ymax) then - surf % i_periodic = i_ymin - else - call fatal_error("Periodic boundary condition applied to & - &interior surface.") - end if - else - surf % i_periodic = surface_dict % get(surf % i_periodic) - end if - - type is (SurfaceZPlane) - if (surf % i_periodic == NONE) then - if (i == i_zmin) then - surf % i_periodic = i_zmax - elseif (i == i_zmax) then - surf % i_periodic = i_zmin - else - call fatal_error("Periodic boundary condition applied to & - &interior surface.") - end if - else - surf % i_periodic = surface_dict % get(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(surf % i_periodic) - end if - - class default - call fatal_error("Periodic boundary condition applied to & - &non-planar surface.") - end select - - ! Make sure opposite surface is also periodic - associate (surf => surfaces(i) % obj) - if (surfaces(surf % i_periodic) % obj % bc /= BC_PERIODIC) then - call fatal_error("Could not find matching surface for periodic & - &boundary on surface " // trim(to_str(surf % id)) // ".") - end if - end associate - end if - end do - ! ========================================================================== ! READ CELLS FROM GEOMETRY.XML diff --git a/src/output.F90 b/src/output.F90 index 1b7cae121e..4113eea472 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)%obj%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 77bded9837..9b82a03967 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -133,12 +133,11 @@ contains real(8), allocatable :: cell_temperatures(:) integer(HID_T) :: geom_group integer(HID_T) :: cells_group, cell_group - integer(HID_T) :: surfaces_group, surface_group + integer(HID_T) :: surfaces_group integer(HID_T) :: universes_group, univ_group integer(HID_T) :: lattices_group, lattice_group character(:), allocatable :: region_spec type(Cell), pointer :: c - class(Surface), pointer :: s class(Lattice), pointer :: lat ! Use H5LT interface to write number of geometry objects @@ -233,7 +232,7 @@ contains region_spec = trim(region_spec) // " |" case default region_spec = trim(region_spec) // " " // to_str(& - sign(surfaces(abs(k))%obj%id, k)) + sign(surfaces(abs(k))%id, k)) end select end do call write_dataset(cell_group, "region", adjustl(region_spec)) diff --git a/src/surface_header.C b/src/surface_header.C index 565352623a..afea2b43d3 100644 --- a/src/surface_header.C +++ b/src/surface_header.C @@ -1,7 +1,8 @@ #include // for std::transform -#include // For std::array +#include #include // For strcmp #include // For numeric_limits +#include #include // For fabs #include "pugixml/pugixml.hpp" @@ -45,6 +46,7 @@ get_node_value(const pugi::xml_node &node, const char *name) std::string value(value_char); std::transform(value.begin(), value.end(), value.begin(), ::tolower); + //TODO: trim whitespace return value; } @@ -57,6 +59,7 @@ get_node_value(const pugi::xml_node &node, const char *name) extern "C" double FP_COINCIDENT; //extern "C" double INFTY; const double INFTY{std::numeric_limits::max()}; +const int C_NONE{-1}; extern "C" const int BC_TRANSMIT{0}; extern "C" const int BC_VACUUM{1}; @@ -70,6 +73,8 @@ extern "C" const int BC_PERIODIC{3}; class Surface; Surface **surfaces_c; +std::map surface_dict; + //============================================================================== // Helper functions for reading the "coeffs" node of an XML surface element //============================================================================== @@ -182,6 +187,19 @@ void read_coeffs(pugi::xml_node surf_node, int surf_id, double &c1, double &c2, } } +//============================================================================== +//============================================================================== + +struct BoundingBox +{ + double xmin; + double xmax; + double ymin; + double ymax; + double zmin; + double zmax; +}; + //============================================================================== //! A geometry primitive used to define regions of 3D space. //============================================================================== @@ -193,8 +211,6 @@ public: //int neighbor_pos[], //!< List of cells on positive side // neighbor_neg[]; //!< List of cells on negative side int bc; //!< Boundary condition - //TODO: switch that zero to a NONE constant. - //int i_periodic = 0; //!< Index of corresponding periodic surface std::string name{""}; //!< User-defined name Surface(pugi::xml_node surf_node); @@ -249,7 +265,6 @@ Surface::Surface(pugi::xml_node surf_node) } else { fatal_error("Must specify id of surface in geometry XML file."); } - //TODO: check for duplicate IDs if (check_for_node(surf_node, "name")) { name = get_node_value(surf_node, "name"); @@ -362,7 +377,9 @@ Surface::to_hdf5(hid_t group_id) const class PeriodicSurface : public Surface { public: - PeriodicSurface(pugi::xml_node surf_node) : Surface(surf_node) {} + int i_periodic{C_NONE}; //!< Index of corresponding periodic surface + + PeriodicSurface(pugi::xml_node surf_node); //! Translate a particle onto this surface from a periodic partner surface. //! @param other A pointer to the partner surface in this periodic BC. @@ -374,8 +391,18 @@ public: //! boundary condition. virtual bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const = 0; + + virtual struct BoundingBox bounding_box() const = 0; }; +PeriodicSurface::PeriodicSurface(pugi::xml_node surf_node) + : Surface(surf_node) +{ + if (check_for_node(surf_node, "periodic_surface_id")) { + i_periodic = stoi(get_node_value(surf_node, "periodic_surface_id")); + } +} + //============================================================================== // Generic functions for x-, y-, and z-, planes. //============================================================================== @@ -428,6 +455,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const; + struct BoundingBox bounding_box() const; }; SurfaceXPlane::SurfaceXPlane(pugi::xml_node surf_node) @@ -485,6 +513,13 @@ bool SurfaceXPlane::periodic_translate(PeriodicSurface *other, double xyz[3], } } +struct BoundingBox +SurfaceXPlane::bounding_box() const +{ + struct BoundingBox out{x0, x0, -INFTY, INFTY, -INFTY, INFTY}; + return out; +} + //============================================================================== // SurfaceYPlane //! A plane perpendicular to the y-axis. @@ -504,6 +539,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const; + struct BoundingBox bounding_box() const; }; SurfaceYPlane::SurfaceYPlane(pugi::xml_node surf_node) @@ -561,6 +597,13 @@ bool SurfaceYPlane::periodic_translate(PeriodicSurface *other, double xyz[3], } } +struct BoundingBox +SurfaceYPlane::bounding_box() const +{ + struct BoundingBox out{-INFTY, INFTY, y0, y0, -INFTY, INFTY}; + return out; +} + //============================================================================== // SurfaceZPlane //! A plane perpendicular to the z-axis. @@ -580,6 +623,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const; + struct BoundingBox bounding_box() const; }; SurfaceZPlane::SurfaceZPlane(pugi::xml_node surf_node) @@ -619,6 +663,13 @@ bool SurfaceZPlane::periodic_translate(PeriodicSurface *other, double xyz[3], return false; } +struct BoundingBox +SurfaceZPlane::bounding_box() const +{ + struct BoundingBox out{-INFTY, INFTY, -INFTY, INFTY, z0, z0}; + return out; +} + //============================================================================== // SurfacePlane //! A general plane. @@ -638,6 +689,7 @@ public: void to_hdf5_inner(hid_t group_id) const; bool periodic_translate(PeriodicSurface *other, double xyz[3], double uvw[3]) const; + struct BoundingBox bounding_box() const; }; SurfacePlane::SurfacePlane(pugi::xml_node surf_node) @@ -698,6 +750,13 @@ bool SurfacePlane::periodic_translate(PeriodicSurface *other, double xyz[3], return false; } +struct BoundingBox +SurfacePlane::bounding_box() const +{ + struct BoundingBox out{-INFTY, INFTY, -INFTY, INFTY, -INFTY, INFTY}; + return out; +} + //============================================================================== // Generic functions for x-, y-, and z-, cylinders //============================================================================== @@ -1359,10 +1418,10 @@ extern "C" void read_surfaces(pugi::xml_node *node) { // Count the number of surfaces. - int n_surfaces = 0; + int n_surfaces{0}; for (pugi::xml_node surf_node = node->child("surface"); surf_node; surf_node = surf_node.next_sibling("surface")) { - n_surfaces += 1; + n_surfaces++; } // Allocate the array of Surface pointers. @@ -1413,12 +1472,125 @@ read_surfaces(pugi::xml_node *node) surfaces_c[i_surf] = new SurfaceQuadric(surf_node); } else { - std::cout << "Call error or handle uppercase here!" << std::endl; + std::cout << "Call error here!" << std::endl; std::cout << surf_type << std::endl; //TODO: call fatal_error } } } + + // Fill the surface dictionary. + for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { + int id = surfaces_c[i_surf]->id; + auto in_dict = surface_dict.find(id); + if (in_dict == surface_dict.end()) { + surface_dict[id] = i_surf; + } else { + std::string err_msg{"Two or more surfaces use the same unique ID: "}; + err_msg += std::to_string(id); + fatal_error(err_msg); + } + } + + // Find the global bounding box (of periodic BC surfaces). + double xmin{INFTY}, xmax{-INFTY}, ymin{INFTY}, ymax{-INFTY}, + zmin{INFTY}, zmax{-INFTY}; + int i_xmin, i_xmax, i_ymin, i_ymax, i_zmin, i_zmax; + for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { + if (surfaces_c[i_surf]->bc == BC_PERIODIC) { + // Downcast to the PeriodicSurface type. + Surface *surf_base = surfaces_c[i_surf]; + PeriodicSurface *surf = dynamic_cast(surf_base); + //TODO: check for null pointer + + // See if this surface makes part of the global bounding box. + struct BoundingBox bb = surf->bounding_box(); + if (bb.xmin > -INFTY and bb.xmin < xmin) { + xmin = bb.xmin; + i_xmin = i_surf; + } + if (bb.xmax < INFTY and bb.xmax > xmax) { + xmax = bb.xmax; + i_xmax = i_surf; + } + if (bb.ymin > -INFTY and bb.ymin < ymin) { + ymin = bb.ymin; + i_ymin = i_surf; + } + if (bb.ymax < INFTY and bb.ymax > ymax) { + ymax = bb.ymax; + i_ymax = i_surf; + } + if (bb.zmin > -INFTY and bb.zmin < zmin) { + zmin = bb.zmin; + i_zmin = i_surf; + } + if (bb.zmax < INFTY and bb.zmax > zmax) { + zmax = bb.zmax; + i_zmax = i_surf; + } + } + } + + // Set i_periodic for periodic BC surfaces. + for (int i_surf = 0; i_surf < n_surfaces; i_surf++) { + if (surfaces_c[i_surf]->bc == BC_PERIODIC) { + // Downcast to the PeriodicSurface type. + Surface *surf_base = surfaces_c[i_surf]; + PeriodicSurface *surf = dynamic_cast(surf_base); + + // Also try downcasting to the SurfacePlane type (which must be handled + // differently). + SurfacePlane *surf_p = dynamic_cast(surf); + + if (surf_p == NULL) { + // This is not a SurfacePlane. + if (surf->i_periodic == C_NONE) { + // The user did not specify the matching periodic surface. See if we + // can find the partnered surface from the bounding box information. + if (i_surf == i_xmin) { + surf->i_periodic = i_xmax; + } else if (i_surf == i_xmax) { + surf->i_periodic = i_xmin; + } else if (i_surf == i_ymin) { + surf->i_periodic = i_ymax; + } else if (i_surf == i_ymax) { + surf->i_periodic = i_ymin; + } else if (i_surf == i_zmin) { + surf->i_periodic = i_zmax; + } else if (i_surf == i_zmax) { + surf->i_periodic = i_zmin; + } else { + fatal_error("Periodic boundary condition applied to interior " + "surface"); + } + } else { + // Convert the surface id to an index. + surf->i_periodic = surface_dict[surf->i_periodic]; + } + } else { + // This is a SurfacePlane. We won't try to find it's partner if the + // user didn't specify one. + if (surf->i_periodic == C_NONE) { + std::string err_msg{"No matching periodic surface specified for " + "periodic boundary condition on surface "}; + err_msg += std::to_string(surf->id); + fatal_error(err_msg); + } else { + // Convert the surface id to an index. + surf->i_periodic = surface_dict[surf->i_periodic]; + } + } + + // Make sure the opposite surface is also periodic. + if (surfaces_c[surf->i_periodic]->bc != BC_PERIODIC) { + std::string err_msg{"Could not find matching surface for periodic " + "boundary condition on surface "}; + err_msg += std::to_string(surf->id); + fatal_error(err_msg); + } + } + } } //============================================================================== @@ -1454,17 +1626,24 @@ surface_to_hdf5(int surf_ind, hid_t group) } extern "C" bool -surface_periodic(int surf_ind1, int surf_ind2, double xyz[3], double uvw[3]) +surface_periodic(int surf_ind1, double xyz[3], double uvw[3]) { - // Hopefully this function has only been called on a pair of surfaces that + // 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[surf_ind2]; + Surface *surf2_gen = surfaces_c[surf1->i_periodic]; PeriodicSurface *surf2 = dynamic_cast(surf2_gen); // Call the type-bound methods. return surf2->periodic_translate(surf1, xyz, uvw); } + +extern "C" int +surface_i_periodic(int surf_ind) +{ + PeriodicSurface *surf = dynamic_cast(surfaces_c[surf_ind]); + return surf->i_periodic; +} diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 98e2423b62..ffe0a6db41 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -12,68 +12,17 @@ module surface_header ! construct closed volumes (cells) !=============================================================================== - type, abstract :: Surface + 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 - integer :: i_periodic = NONE ! Index of corresponding periodic surface - character(len=104) :: name = "" ! User-defined name end type Surface -!=============================================================================== -! SURFACECONTAINER allows us to store an array of different types of surfaces -!=============================================================================== - - type :: SurfaceContainer - class(Surface), allocatable :: obj - end type SurfaceContainer - -!=============================================================================== -! All the derived types below are extensions of the abstract Surface type. They -! inherent the reflect() type-bound procedure and must implement normal() -!=============================================================================== - - type, extends(Surface) :: SurfaceXPlane - end type SurfaceXPlane - - type, extends(Surface) :: SurfaceYPlane - end type SurfaceYPlane - - type, extends(Surface) :: SurfaceZPlane - end type SurfaceZPlane - - type, extends(Surface) :: SurfacePlane - end type SurfacePlane - - type, extends(Surface) :: SurfaceXCylinder - end type SurfaceXCylinder - - type, extends(Surface) :: SurfaceYCylinder - end type SurfaceYCylinder - - type, extends(Surface) :: SurfaceZCylinder - end type SurfaceZCylinder - - type, extends(Surface) :: SurfaceSphere - end type SurfaceSphere - - type, extends(Surface) :: SurfaceXCone - end type SurfaceXCone - - type, extends(Surface) :: SurfaceYCone - end type SurfaceYCone - - type, extends(Surface) :: SurfaceZCone - end type SurfaceZCone - - type, extends(Surface) :: SurfaceQuadric - end type SurfaceQuadric - integer(C_INT32_T), bind(C) :: n_surfaces ! # of surfaces - type(SurfaceContainer), allocatable, target :: surfaces(:) + type(Surface), allocatable, target :: surfaces(:) ! Dictionary that maps user IDs to indices in 'surfaces' type(DictIntInt) :: surface_dict diff --git a/src/tallies/tally_filter_surface.F90 b/src/tallies/tally_filter_surface.F90 index daecd88917..6a2afa3ec6 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)) % obj % 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 d616212ec6..98bca2a279 100644 --- a/src/tracking.F90 +++ b/src/tracking.F90 @@ -6,7 +6,7 @@ module tracking 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_periodic_c, surface_i_periodic_c use message_passing use mgxs_header use nuclide_header @@ -298,7 +298,7 @@ contains class(Surface), pointer :: surf i_surface = abs(p % surface) - surf => surfaces(i_surface)%obj + surf => surfaces(i_surface) if (verbosity >= 10 .or. trace) then call write_message(" Crossing surface " // trim(to_str(surf % id))) end if @@ -412,14 +412,14 @@ contains p % coord(1) % xyz = xyz end if - rotational = surface_periodic_c(i_surface-1, surf % i_periodic-1, & - p % coord(1) % xyz, p % coord(1) % uvw) + rotational = surface_periodic_c(i_surface-1, p % coord(1) % xyz, & + p % coord(1) % uvw) ! Reassign particle's surface if (rotational) then - p % surface = surf % i_periodic + p % surface = surface_i_periodic_c(i_surface-1) + 1 else - p % surface = sign(surf % i_periodic, p % surface) + p % surface = sign(surface_i_periodic_c(i_surface-1) + 1, p % surface) end if ! Figure out what cell particle is in now