diff --git a/src/cell.cpp b/src/cell.cpp index c6f7b3ab9e..83516eba67 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -748,8 +748,6 @@ extern "C" { int32_t cell_offset(Cell* c, int map) {return c->offset[map];} - void cell_to_hdf5(Cell* c, hid_t group) {c->to_hdf5(group);} - void extend_cells_c(int32_t n) { global_cells.reserve(global_cells.size() + n); diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index df59e34101..b58b5415ba 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -111,12 +111,6 @@ module geometry_header integer(C_INT32_T) :: offset end function cell_offset_c - subroutine cell_to_hdf5_c(cell_ptr, group) bind(C, name='cell_to_hdf5') - import HID_T, C_PTR - type(C_PTR), intent(in), value :: cell_ptr - integer(HID_T), intent(in), value :: group - end subroutine cell_to_hdf5_c - function lattice_pointer(lat_ind) bind(C) result(ptr) import C_PTR, C_INT32_T integer(C_INT32_T), intent(in), value :: lat_ind @@ -137,29 +131,6 @@ module geometry_header logical(C_BOOL) :: is_valid end function lattice_are_valid_indices_c - subroutine lattice_get_indices_c(lat_ptr, xyz, i_xyz) & - bind(C, name='lattice_get_indices') - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), intent(in), value :: lat_ptr - real(C_DOUBLE), intent(in) :: xyz(3) - integer(C_INT), intent(out) :: i_xyz(3) - end subroutine lattice_get_indices_c - - subroutine lattice_get_local_xyz_c(lat_ptr, global_xyz, i_xyz, local_xyz) & - bind(C, name='lattice_get_local_xyz') - import C_PTR, C_INT, C_DOUBLE - type(C_PTR), intent(in), value :: lat_ptr - real(C_DOUBLE), intent(in) :: global_xyz(3) - integer(C_INT), intent(in) :: i_xyz(3) - real(C_DOUBLE), intent(out) :: local_xyz(3) - end subroutine lattice_get_local_xyz_c - - subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5') - import HID_T, C_PTR - type(C_PTR), intent(in), value :: lat_ptr - integer(HID_T), intent(in), value :: group - end subroutine lattice_to_hdf5_c - function lattice_offset_c(lat_ptr, map, i_xyz) & bind(C, name='lattice_offset') result(offset) import C_PTR, C_INT, C_INT32_T @@ -169,14 +140,6 @@ module geometry_header integer(C_INT32_T) :: offset end function lattice_offset_c - function lattice_universe_c(lat_ptr, i_xyz) & - bind(C, name='lattice_universe') result(univ) - import C_PTR, C_INT32_t, C_INT - type(C_PTR), intent(in), value :: lat_ptr - integer(C_INT), intent(in) :: i_xyz(3) - integer(C_INT32_T) :: univ - end function lattice_universe_c - subroutine extend_cells_c(n) bind(C) import C_INT32_t integer(C_INT32_T), intent(in), value :: n @@ -187,40 +150,14 @@ module geometry_header ! LATTICE abstract type for ordered array of universes. !=============================================================================== - type, abstract :: Lattice + type :: Lattice type(C_PTR) :: ptr contains procedure :: id => lattice_id procedure :: are_valid_indices => lattice_are_valid_indices - procedure :: get => lattice_get - procedure :: get_indices => lattice_get_indices - procedure :: get_local_xyz => lattice_get_local_xyz procedure :: offset => lattice_offset - procedure :: to_hdf5 => lattice_to_hdf5 end type Lattice -!=============================================================================== -! RECTLATTICE extends LATTICE for rectilinear arrays. -!=============================================================================== - - type, extends(Lattice) :: RectLattice - end type RectLattice - -!=============================================================================== -! HEXLATTICE extends LATTICE for hexagonal (sometimes called triangular) arrays. -!=============================================================================== - - type, extends(Lattice) :: HexLattice - end type HexLattice - -!=============================================================================== -! LATTICECONTAINER pointer array for storing lattices -!=============================================================================== - - type LatticeContainer - class(Lattice), allocatable :: obj - end type LatticeContainer - !=============================================================================== ! CELL defines a closed volume by its bounding surfaces !=============================================================================== @@ -246,7 +183,6 @@ module geometry_header procedure :: sqrtkT_size => cell_sqrtkT_size procedure :: sqrtkT => cell_sqrtkT procedure :: offset => cell_offset - procedure :: to_hdf5 => cell_to_hdf5 end type Cell @@ -257,7 +193,7 @@ module geometry_header integer(C_INT32_T), bind(C) :: n_universes ! # of universes type(Cell), allocatable, target :: cells(:) - type(LatticeContainer), allocatable, target :: lattices(:) + type(Lattice), allocatable, target :: lattices(:) ! Dictionaries which map user IDs to indices in the global arrays type(DictIntInt) :: cell_dict @@ -279,29 +215,6 @@ contains is_valid = lattice_are_valid_indices_c(this % ptr, i_xyz) end function lattice_are_valid_indices - function lattice_get(this, i_xyz) result(univ) - class(Lattice), intent(in) :: this - integer(C_INT), intent(in) :: i_xyz(3) - integer(C_INT32_T) :: univ - univ = lattice_universe_c(this % ptr, i_xyz) - end function lattice_get - - function lattice_get_indices(this, xyz) result(i_xyz) - class(Lattice), intent(in) :: this - real(C_DOUBLE), intent(in) :: xyz(3) - integer(C_INT) :: i_xyz(3) - call lattice_get_indices_c(this % ptr, xyz, i_xyz) - end function lattice_get_indices - - function lattice_get_local_xyz(this, global_xyz, i_xyz) & - result(local_xyz) - class(Lattice), intent(in) :: this - real(C_DOUBLE), intent(in) :: global_xyz(3) - integer(C_INT), intent(in) :: i_xyz(3) - real(C_DOUBLE) :: local_xyz(3) - call lattice_get_local_xyz_c(this % ptr, global_xyz, i_xyz, local_xyz) - end function lattice_get_local_xyz - function lattice_offset(this, map, i_xyz) result(offset) class(Lattice), intent(in) :: this integer(C_INT), intent(in) :: map @@ -310,12 +223,6 @@ contains offset = lattice_offset_c(this % ptr, map, i_xyz) end function lattice_offset - subroutine lattice_to_hdf5(this, group) - class(Lattice), intent(in) :: this - integer(HID_T), intent(in) :: group - call lattice_to_hdf5_c(this % ptr, group) - end subroutine lattice_to_hdf5 - !=============================================================================== function cell_id(this) result(id) @@ -393,12 +300,6 @@ contains offset = cell_offset_c(this % ptr, map) end function cell_offset - subroutine cell_to_hdf5(this, group) - class(Cell), intent(in) :: this - integer(HID_T), intent(in) :: group - call cell_to_hdf5_c(this % ptr, group) - end subroutine cell_to_hdf5 - !=============================================================================== ! GET_TEMPERATURES returns a list of temperatures that each nuclide/S(a,b) table ! appears at in the model. Later, this list is used to determine the actual diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 1a9a234f58..cbdc54fec9 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1170,11 +1170,8 @@ contains allocate(lattices(n_rlats + n_hlats)) RECT_LATTICES: do i = 1, n_rlats - allocate(RectLattice::lattices(i) % obj) - lat => lattices(i) % obj + lat => lattices(i) lat % ptr = lattice_pointer(i - 1) - select type(lat) - type is (RectLattice) ! Get pointer to i-th lattice node_lat = node_rlat_list(i) @@ -1182,15 +1179,11 @@ contains ! Add lattice to dictionary call lattice_dict % set(lat % id(), i) - end select end do RECT_LATTICES HEX_LATTICES: do i = 1, n_hlats - allocate(HexLattice::lattices(n_rlats + i) % obj) - lat => lattices(n_rlats + i) % obj + lat => lattices(n_rlats + i) lat % ptr = lattice_pointer(n_rlats + i - 1) - select type (lat) - type is (HexLattice) ! Get pointer to i-th lattice node_lat = node_hlat_list(i) @@ -1198,7 +1191,6 @@ contains ! Add lattice to dictionary call lattice_dict % set(lat % id(), n_rlats + i) - end select end do HEX_LATTICES ! ========================================================================== diff --git a/src/lattice.cpp b/src/lattice.cpp index 46f88e3519..241c4dedcf 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -206,7 +206,7 @@ RectLattice::RectLattice(pugi::xml_node lat_node) //============================================================================== int32_t& -RectLattice::operator[](const int i_xyz[3]) +RectLattice::operator[](std::array i_xyz) { int indx = nx*ny*i_xyz[2] + nx*i_xyz[1] + i_xyz[0]; return universes[indx]; @@ -225,7 +225,8 @@ RectLattice::are_valid_indices(const int i_xyz[3]) const //============================================================================== std::pair> -RectLattice::distance(Position r, Direction u, const int i_xyz[3]) const +RectLattice::distance(Position r, Direction u, const std::array& i_xyz) +const { // Get short aliases to the coordinates. double x = r.x; @@ -299,7 +300,8 @@ RectLattice::get_indices(Position r) const //============================================================================== Position -RectLattice::get_local_position(Position r, const int i_xyz[3]) const +RectLattice::get_local_position(Position r, const std::array i_xyz) +const { r.x -= (lower_left.x + (i_xyz[0] + 0.5)*pitch.x); r.y -= (lower_left.y + (i_xyz[1] + 0.5)*pitch.y); @@ -549,7 +551,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node) //============================================================================== int32_t& -HexLattice::operator[](const int i_xyz[3]) +HexLattice::operator[](std::array i_xyz) { int indx = (2*n_rings-1)*(2*n_rings-1) * i_xyz[2] + (2*n_rings-1) * i_xyz[1] @@ -580,7 +582,8 @@ HexLattice::are_valid_indices(const int i_xyz[3]) const //============================================================================== std::pair> -HexLattice::distance(Position r, Direction u, const int i_xyz[3]) const +HexLattice::distance(Position r, Direction u, const std::array& i_xyz) +const { // Compute the direction on the hexagonal basis. double beta_dir = u.x * std::sqrt(3.0) / 2.0 + u.y / 2.0; @@ -598,10 +601,10 @@ HexLattice::distance(Position r, Direction u, const int i_xyz[3]) const double edge = -copysign(0.5*pitch[0], beta_dir); // Oncoming edge Position r_t; if (beta_dir > 0) { - const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; + const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; r_t = get_local_position(r, i_xyz_t); } else { - const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; + const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; r_t = get_local_position(r, i_xyz_t); } double beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0; @@ -617,10 +620,10 @@ HexLattice::distance(Position r, Direction u, const int i_xyz[3]) const // Lower-right and upper-left sides. edge = -copysign(0.5*pitch[0], gamma_dir); if (gamma_dir > 0) { - const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; + const std::array i_xyz_t {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; r_t = get_local_position(r, i_xyz_t); } else { - const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; + const std::array i_xyz_t {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; r_t = get_local_position(r, i_xyz_t); } double gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0; @@ -639,10 +642,10 @@ HexLattice::distance(Position r, Direction u, const int i_xyz[3]) const // Upper and lower sides. edge = -copysign(0.5*pitch[0], u.y); if (u.y > 0) { - const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; + const std::array i_xyz_t {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; r_t = get_local_position(r, i_xyz_t); } else { - const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; + const std::array i_xyz_t {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; r_t = get_local_position(r, i_xyz_t); } if ((std::abs(r_t.y - edge) > FP_PRECISION) && u.y != 0) { @@ -719,7 +722,7 @@ HexLattice::get_indices(Position r) const double d_min {INFTY}; for (int i = 0; i < 2; i++) { for (int j = 0; j < 2; j++) { - int i_xyz[3] {out[0] + j, out[1] + i, 0}; + const std::array i_xyz {out[0] + j, out[1] + i, 0}; Position r_t = get_local_position(r, i_xyz); double d = r_t.x*r_t.x + r_t.y*r_t.y; if (d < d_min) { @@ -747,7 +750,8 @@ HexLattice::get_indices(Position r) const //============================================================================== Position -HexLattice::get_local_position(Position r, const int i_xyz[3]) const +HexLattice::get_local_position(Position r, const std::array i_xyz) +const { // x_l = x_g - (center + pitch_x*cos(30)*index_x) r.x -= (center.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]); @@ -893,32 +897,8 @@ extern "C" { bool lattice_are_valid_indices(Lattice *lat, const int i_xyz[3]) {return lat->are_valid_indices(i_xyz);} - void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3]) - { - Position r {xyz}; - std::array inds = lat->get_indices(r); - i_xyz[0] = inds[0]; - i_xyz[1] = inds[1]; - i_xyz[2] = inds[2]; - } - - void lattice_get_local_xyz(Lattice *lat, const double global_xyz[3], - const int i_xyz[3], double local_xyz[3]) - { - Position global {global_xyz}; - Position local = lat->get_local_position(global, i_xyz); - local_xyz[0] = local.x; - local_xyz[1] = local.y; - local_xyz[2] = local.z; - } - int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3]) {return lat->offset(map, i_xyz);} - - void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);} - - int32_t lattice_universe(Lattice *lat, const int i_xyz[3]) - {return (*lat)[i_xyz];} } } // namespace openmc diff --git a/src/lattice.h b/src/lattice.h index 758acfa829..163b86e851 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -57,14 +57,7 @@ public: virtual ~Lattice() {} - virtual int32_t& operator[](const int i_xyz[3]) = 0; - - int32_t& - operator[](std::array i_xyz) - { - int i_xyz_[3] {i_xyz[0], i_xyz[1], i_xyz[2]}; - return operator[](i_xyz_); - } + virtual int32_t& operator[](std::array i_xyz) = 0; virtual LatticeIter begin(); LatticeIter end(); @@ -98,20 +91,13 @@ public: //! \brief Find the next lattice surface crossing //! \param r A 3D Cartesian coordinate. //! \param u A 3D Cartesian direction. - //! \param i_xyz[3] The indices for a lattice tile. + //! \param i_xyz The indices for a lattice tile. //! \return The distance to the next crossing and an array indicating how the //! lattice indices would change after crossing that boundary. virtual std::pair> - distance(Position r, Direction u, const int i_xyz[3]) const + distance(Position r, Direction u, const std::array& i_xyz) const = 0; - std::pair> - distance(Position r, Direction u, std::array i_xyz) const - { - int i_xyz_[3] {i_xyz[0], i_xyz[1], i_xyz[2]}; - return distance(r, u, i_xyz_); - } - //! \brief Find the lattice tile indices for a given point. //! \param r A 3D Cartesian coordinate. //! \return An array containing the indices of a lattice tile. @@ -119,17 +105,10 @@ public: //! \brief Get coordinates local to a lattice tile. //! \param r A 3D Cartesian coordinate. - //! \param i_xyz[3] The indices for a lattice tile. + //! \param i_xyz The indices for a lattice tile. //! \return Local 3D Cartesian coordinates. virtual Position - get_local_position(Position r, const int i_xyz[3]) const = 0; - - Position - get_local_position(Position r, std::array i_xyz) const - { - int i_xyz_[3] {i_xyz[0], i_xyz[1], i_xyz[2]}; - return get_local_position(r, i_xyz_); - } + get_local_position(Position r, const std::array i_xyz) const = 0; //! \brief Check flattened lattice index. //! \param indx The index for a lattice tile. @@ -223,17 +202,17 @@ class RectLattice : public Lattice public: explicit RectLattice(pugi::xml_node lat_node); - int32_t& operator[](const int i_xyz[3]); + int32_t& operator[](std::array i_xyz); bool are_valid_indices(const int i_xyz[3]) const; std::pair> - distance(Position r, Direction u, const int i_xyz[3]) const; + distance(Position r, Direction u, const std::array& i_xyz) const; std::array get_indices(Position r) const; Position - get_local_position(Position r, const int i_xyz[3]) const; + get_local_position(Position r, const std::array i_xyz) const; int32_t& offset(int map, const int i_xyz[3]); @@ -259,7 +238,7 @@ class HexLattice : public Lattice public: explicit HexLattice(pugi::xml_node lat_node); - int32_t& operator[](const int i_xyz[3]); + int32_t& operator[](std::array i_xyz); LatticeIter begin(); @@ -268,12 +247,12 @@ public: bool are_valid_indices(const int i_xyz[3]) const; std::pair> - distance(Position r, Direction u, const int i_xyz[3]) const; + distance(Position r, Direction u, const std::array& i_xyz) const; std::array get_indices(Position r) const; Position - get_local_position(Position r, const int i_xyz[3]) const; + get_local_position(Position r, const std::array i_xyz) const; bool is_valid_index(int indx) const; diff --git a/src/output.F90 b/src/output.F90 index 0f23c67bd8..a4d7341542 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -244,7 +244,7 @@ contains ! Print information on lattice if (p % coord(i) % lattice /= NONE) then - l => lattices(p % coord(i) % lattice) % obj + l => lattices(p % coord(i) % lattice) write(ou,*) ' Lattice = ' // trim(to_str(l % id())) write(ou,*) ' Lattice position = (' // trim(to_str(& p % coord(i) % lattice_x)) // ',' // trim(to_str(& diff --git a/src/surface.cpp b/src/surface.cpp index 421c43d485..5ae3900357 100644 --- a/src/surface.cpp +++ b/src/surface.cpp @@ -1228,17 +1228,6 @@ extern "C" { uvw[2] = u.z; } - void surface_normal(Surface* surf, double xyz[3], double uvw[3]) - { - Position r {xyz}; - Direction u = surf->normal(r); - uvw[0] = u.x; - uvw[1] = u.y; - uvw[2] = u.z; - } - - void surface_to_hdf5(Surface* surf, hid_t group) {surf->to_hdf5(group);} - int surface_i_periodic(PeriodicSurface* surf) {return surf->i_periodic;} bool diff --git a/src/surface_header.F90 b/src/surface_header.F90 index 401ac90aae..9d7e114613 100644 --- a/src/surface_header.F90 +++ b/src/surface_header.F90 @@ -38,23 +38,6 @@ module surface_header real(C_DOUBLE), intent(inout) :: uvw(3); end subroutine surface_reflect_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') - import C_PTR, HID_T - 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 @@ -91,8 +74,6 @@ module surface_header procedure :: id => surface_id procedure :: bc => surface_bc procedure :: reflect => surface_reflect - procedure :: normal => surface_normal - procedure :: to_hdf5 => surface_to_hdf5 procedure :: i_periodic => surface_i_periodic procedure :: periodic_translate => surface_periodic @@ -126,19 +107,6 @@ contains call surface_reflect_c(this % ptr, xyz, uvw) end subroutine surface_reflect - 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 diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index 48b223a439..ae4e09023d 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -61,12 +61,11 @@ contains offset = offset + cells(p % coord(i) % cell + 1) & % offset(distribcell_index-1) elseif (cells(p % coord(i) % cell + 1) % type() == FILL_LATTICE) then - if (lattices(p % coord(i + 1) % lattice) % obj & - % are_valid_indices([& + if (lattices(p % coord(i + 1) % lattice) % are_valid_indices([& p % coord(i + 1) % lattice_x, & p % coord(i + 1) % lattice_y, & p % coord(i + 1) % lattice_z])) then - offset = offset + lattices(p % coord(i + 1) % lattice) % obj & + offset = offset + lattices(p % coord(i + 1) % lattice) & % offset(distribcell_index - 1, & [p % coord(i + 1) % lattice_x, & p % coord(i + 1) % lattice_y, &