diff --git a/src/geometry.F90 b/src/geometry.F90 index f717810bb1..29e03bbb5f 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -69,9 +69,7 @@ contains type(Particle), intent(inout) :: p logical, intent(inout) :: found integer, optional :: search_cells(:) - integer :: j, k ! coordinate level index - integer :: offset ! instance # of a distributed cell - integer :: distribcell_index + integer :: j ! coordinate level index integer :: i_xyz(3) ! indices in lattice integer :: i_cell ! index in cells array @@ -85,61 +83,7 @@ contains if (found) then associate(c => cells(i_cell)) - CELL_TYPE: if (c % type() == FILL_MATERIAL) then - ! ====================================================================== - ! AT LOWEST UNIVERSE, TERMINATE SEARCH - - ! Save previous material and temperature - p % last_material = p % material - p % last_sqrtkT = p % sqrtkT - - ! Get distributed offset - if (c % material_size() > 1 .or. c % sqrtkT_size() > 1) then - ! Distributed instances of this cell have different - ! materials/temperatures. Determine which instance this is for - ! assigning the matching material/temperature. - distribcell_index = c % distribcell_index() - offset = 0 - do k = 1, p % n_coord - if (cells(p % coord(k) % cell) % type() == FILL_UNIVERSE) then - offset = offset + cells(p % coord(k) % cell) & - % offset(distribcell_index-1) - elseif (cells(p % coord(k) % cell) % type() == FILL_LATTICE) then - if (lattices(p % coord(k + 1) % lattice) % obj & - % are_valid_indices([& - p % coord(k + 1) % lattice_x, & - p % coord(k + 1) % lattice_y, & - p % coord(k + 1) % lattice_z])) then - offset = offset + lattices(p % coord(k + 1) % lattice) % obj & - % offset(distribcell_index - 1, & - [p % coord(k + 1) % lattice_x, & - p % coord(k + 1) % lattice_y, & - p % coord(k + 1) % lattice_z]) - end if - end if - end do - - ! Keep track of which instance of the cell the particle is in - p % cell_instance = offset + 1 - else - p % cell_instance = 1 - end if - - ! Save the material - if (c % material_size() > 1) then - p % material = c % material(offset + 1) - else - p % material = c % material(1) - end if - - ! Save the temperature - if (c % sqrtkT_size() > 1) then - p % sqrtkT = c % sqrtkT(offset) - else - p % sqrtkT = c % sqrtkT(0) - end if - - elseif (c % type() == FILL_UNIVERSE) then CELL_TYPE + CELL_TYPE: if (c % type() == FILL_UNIVERSE) then ! ====================================================================== ! CELL CONTAINS LOWER UNIVERSE, RECURSIVELY FIND CELL diff --git a/src/geometry.cpp b/src/geometry.cpp index 4c31fe85d0..7131a45135 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -5,7 +5,7 @@ #include "cell.h" #include "constants.h" #include "error.h" -#include "particle.h" +#include "lattice.h" #include "settings.h" //TODO: remove this include @@ -16,6 +16,8 @@ namespace openmc { std::vector overlap_check_count; +//============================================================================== + extern "C" bool check_cell_overlap(Particle* p) { int n_coord = p->n_coord; @@ -72,10 +74,11 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) { // Find which cell of this universe the particle is in. bool found = false; + int32_t i_cell; for (int i = 0; i < n_search_cells; i++) { - int32_t i_cell = search_cells[i]; + i_cell = search_cells[i]; - // Make sure the search cell is in the same universe + // Make sure the search cell is in the same universe. //TODO: off-by-one indexing if (global_cells[i_cell]->universe - 1 != i_universe) continue; @@ -96,6 +99,58 @@ find_cell(Particle* p, int n_search_cells, int* search_cells) { } } + if (found) { + Cell& c {*global_cells[i_cell]}; + if (c.type == FILL_MATERIAL) { + // Find the distribcell instance number. + if (c.material.size() > 1 || c.sqrtkT.size() > 1) { + //===================================================================== + //! Found a material cell which means this is the lowest coord level. + + //TODO: off-by-one indexing + int distribcell_index = c.distribcell_index - 1; + int offset = 0; + for (int i = 0; i < p->n_coord; i++) { + Cell& c_i {*global_cells[p->coord[i].cell-1]}; + if (c_i.type == FILL_UNIVERSE) { + offset += c_i.offset[distribcell_index]; + } else if (c_i.type == FILL_LATTICE) { + Lattice& lat {*lattices_c[p->coord[i+1].lattice-1]}; + int i_xyz[3] {p->coord[i+1].lattice_x, + p->coord[i+1].lattice_y, + p->coord[i+1].lattice_z}; + if (lat.are_valid_indices(i_xyz)) { + offset += lat.offset(distribcell_index, i_xyz); + } + } + } + p->cell_instance = offset + 1; + } else { + p->cell_instance = 1; + } + + // Set the material and temperature. + p->last_material = p->material; + int32_t mat; + if (c.material.size() > 1) { + mat = c.material[p->cell_instance-1]; + } else { + mat = c.material[0]; + } + if (mat == MATERIAL_VOID) { + p->material = MATERIAL_VOID; + } else { + p->material = mat + 1; + } + p->last_sqrtkT = p->sqrtkT; + if (c.sqrtkT.size() > 1) { + p->sqrtkT = c.sqrtkT[p->cell_instance-1]; + } else { + p->sqrtkT = c.sqrtkT[0]; + } + } + } + return found; } diff --git a/src/geometry.h b/src/geometry.h index e02a17b983..b73ae314c4 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -4,13 +4,28 @@ #include #include +#include "particle.h" + namespace openmc { extern "C" int openmc_root_universe; -//TODO: free this memory extern std::vector overlap_check_count; +//============================================================================== +//! Check for overlapping cells at the particle's position. +//============================================================================== + +extern "C" bool +check_cell_overlap(Particle* p); + +//============================================================================== +//! Locate the particle in the geometry tree and set its geometry data fields. +//============================================================================== + +extern "C" bool +find_cell(Particle* p, int n_search_cells, int* search_cells); + } // namespace openmc #endif // OPENMC_GEOMETRY_H