diff --git a/src/geometry.F90 b/src/geometry.F90 index 817da18de..dcaabb9dd 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -65,78 +65,15 @@ contains !=============================================================================== recursive subroutine find_cell(p, found, search_cells) - type(Particle), intent(inout) :: p logical, intent(inout) :: found integer, optional :: search_cells(:) - integer :: j ! coordinate level index - integer :: i_xyz(3) ! indices in lattice - integer :: i_cell ! index in cells array if (present(search_cells)) then found = find_cell_c(p, size(search_cells), search_cells-1) else found = find_cell_c(p, 0) end if - j = p % n_coord - i_cell = p % coord(j) % cell - - if (found) then - associate(c => cells(i_cell)) - CELL_TYPE: if (c % type() == FILL_UNIVERSE) then - j = j + 1 - p % n_coord = j - call find_cell(p, found) - j = p % n_coord - - elseif (c % type() == FILL_LATTICE) then CELL_TYPE - ! ====================================================================== - ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - - associate (lat => lattices(c % fill() + 1) % obj) - ! Determine lattice indices - i_xyz = lat % get_indices(p % coord(j) % xyz + TINY_BIT * p % coord(j) % uvw) - - ! Store lower level coordinates - p % coord(j + 1) % xyz = lat % get_local_xyz(p % coord(j) % xyz, i_xyz) - p % coord(j + 1) % uvw = p % coord(j) % uvw - - ! set particle lattice indices - p % coord(j + 1) % lattice = c % fill() + 1 - p % coord(j + 1) % lattice_x = i_xyz(1) - p % coord(j + 1) % lattice_y = i_xyz(2) - p % coord(j + 1) % lattice_z = i_xyz(3) - - ! Set the next lowest coordinate level. - if (lat % are_valid_indices(i_xyz)) then - ! Particle is inside the lattice. - p % coord(j + 1) % universe = & - lat % get([i_xyz(1), i_xyz(2), i_xyz(3)]) + 1 - - else - ! Particle is outside the lattice. - if (lat % outer() == NO_OUTER_UNIVERSE) then - call warning("Particle " // trim(to_str(p %id)) & - // " is outside lattice " // trim(to_str(lat % id())) & - // " but the lattice has no defined outer universe.") - found = .false. - return - else - p % coord(j + 1) % universe = lat % outer() + 1 - end if - end if - end associate - - ! Move particle to next level and search for the lower cells. - j = j + 1 - p % n_coord = j - - call find_cell(p, found) - j = p % n_coord - - end if CELL_TYPE - end associate - end if end subroutine find_cell diff --git a/src/geometry.cpp b/src/geometry.cpp index 695afe0be..e400a2245 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -8,9 +8,6 @@ #include "lattice.h" #include "settings.h" -//TODO: remove this include -#include - namespace openmc { @@ -153,8 +150,7 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) { //======================================================================== //! Found a lower universe, update this coord level then search the next. - // Add another coordinate level. - //++p->n_coord; + // Set the lower coordinate level universe. p->coord[p->n_coord].universe = c.fill + 1; // Set the position and direction. @@ -190,6 +186,57 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) { + w*c.rotation[11]; p->coord[p->n_coord].rotated = true; } + + // Update the coordinate level and recurse. + ++p->n_coord; + find_cell(p, 0, nullptr); + + } else if (c.type == FILL_LATTICE) { + //======================================================================== + //! Found a lower lattice, update this coord level then search the next. + + Lattice& lat {*lattices_c[c.fill]}; + + // Determine lattice indices. + Position r {p->coord[p->n_coord-1].xyz}; + Direction u {p->coord[p->n_coord-1].uvw}; + r += TINY_BIT * u; + auto i_xyz = lat.get_indices(r); + + // Store lower level coordinates. + r = lat.get_local_position(p->coord[p->n_coord-1].xyz, i_xyz); + p->coord[p->n_coord].xyz[0] = r.x; + p->coord[p->n_coord].xyz[1] = r.y; + p->coord[p->n_coord].xyz[2] = r.z; + p->coord[p->n_coord].uvw[0] = u.x; + p->coord[p->n_coord].uvw[1] = u.y; + p->coord[p->n_coord].uvw[2] = u.z; + + // Set lattice indices. + p->coord[p->n_coord].lattice = c.fill + 1; + p->coord[p->n_coord].lattice_x = i_xyz[0]; + p->coord[p->n_coord].lattice_y = i_xyz[1]; + p->coord[p->n_coord].lattice_z = i_xyz[2]; + + // Set the lower coordinate level universe. + if (lat.are_valid_indices(i_xyz)) { + p->coord[p->n_coord].universe = lat[i_xyz] + 1; + } else { + if (lat.outer != NO_OUTER_UNIVERSE) { + p->coord[p->n_coord].universe = lat.outer + 1; + } else { + std::stringstream err_msg; + err_msg << "Particle " << p->id << " is outside lattice " + << lat.id << " but the lattice has no defined outer " + "universe."; + warning(err_msg); + found = false; + } + } + + // Update the coordinate level and recurse. + ++p->n_coord; + find_cell(p, 0, nullptr); } } diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 272966b33..23700d829 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -191,13 +191,6 @@ module geometry_header integer(C_INT32_T) :: offset end function lattice_offset_c - function lattice_outer_c(lat_ptr) bind(C, name='lattice_outer') & - result(outer) - import C_PTR, C_INT32_T - type(C_PTR), intent(in), value :: lat_ptr - integer(C_INT32_T) :: outer - end function lattice_outer_c - function lattice_universe_c(lat_ptr, i_xyz) & bind(C, name='lattice_universe') result(univ) import C_PTR, C_INT32_t, C_INT @@ -238,7 +231,6 @@ module geometry_header procedure :: get_indices => lattice_get_indices procedure :: get_local_xyz => lattice_get_local_xyz procedure :: offset => lattice_offset - procedure :: outer => lattice_outer procedure :: to_hdf5 => lattice_to_hdf5 end type Lattice @@ -275,7 +267,6 @@ module geometry_header ! Boolean expression of half-spaces ! Rotation matrix and translation vector - real(8), allocatable :: translation(:) real(8), allocatable :: rotation(:) real(8), allocatable :: rotation_matrix(:,:) @@ -370,12 +361,6 @@ contains offset = lattice_offset_c(this % ptr, map, i_xyz) end function lattice_offset - function lattice_outer(this) result(outer) - class(Lattice), intent(in) :: this - integer(C_INT32_T) :: outer - outer = lattice_outer_c(this % ptr) - end function lattice_outer - subroutine lattice_to_hdf5(this, group) class(Lattice), intent(in) :: this integer(HID_T), intent(in) :: group diff --git a/src/lattice.cpp b/src/lattice.cpp index e2871d9e1..8797d0fba 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -923,8 +923,6 @@ extern "C" { int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3]) {return lat->offset(map, i_xyz);} - int32_t lattice_outer(Lattice *lat) {return lat->outer;} - void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);} int32_t lattice_universe(Lattice *lat, const int i_xyz[3]) diff --git a/src/lattice.h b/src/lattice.h index 94d7b4387..df7c699d8 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -54,6 +54,13 @@ public: 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 LatticeIter begin(); LatticeIter end(); @@ -76,6 +83,13 @@ public: //! otherwise. virtual bool are_valid_indices(const int i_xyz[3]) const = 0; + bool + are_valid_indices(std::array i_xyz) const + { + int i_xyz_[3] {i_xyz[0], i_xyz[1], i_xyz[2]}; + return are_valid_indices(i_xyz_); + } + //! \brief Find the next lattice surface crossing //! \param r A 3D Cartesian coordinate. //! \param u A 3D Cartesian direction. @@ -98,6 +112,13 @@ public: 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_); + } + //! \brief Check flattened lattice index. //! \param indx The index for a lattice tile. //! \return true if the given index fit within the lattice bounds. False