diff --git a/src/cell.cpp b/src/cell.cpp index 1d630b3d4..82138a83f 100644 --- a/src/cell.cpp +++ b/src/cell.cpp @@ -234,6 +234,21 @@ Cell::Cell(pugi::xml_node cell_node) material.shrink_to_fit(); } + // Make sure that either material or fill was specified. + if ((material[0] == C_NONE) && (fill == C_NONE)) { + std::stringstream err_msg; + err_msg << "Neither material nor fill was specified for cell " << id; + fatal_error(err_msg); + } + + // Make sure that material and fill haven't been specified simultaneously. + if ((material[0] != C_NONE) && (fill != C_NONE)) { + std::stringstream err_msg; + err_msg << "Cell " << id << " has both a material and a fill specified; " + << "only one can be specified per cell"; + fatal_error(err_msg); + } + std::string region_spec {""}; if (check_for_node(cell_node, "region")) { region_spec = get_node_value(cell_node, "region"); @@ -470,6 +485,10 @@ extern "C" { void cell_set_universe(Cell *c, int32_t universe) {c->universe = universe;} + int32_t cell_fill(Cell *c) {return c->fill;} + + int32_t* cell_fill_ptr(Cell *c) {return &c->fill;} + int32_t cell_n_instances(Cell *c) {return c->n_instances;} bool cell_simple(Cell *c) {return c->simple;} diff --git a/src/constants.F90 b/src/constants.F90 index dbd4181da..d75ee4e54 100644 --- a/src/constants.F90 +++ b/src/constants.F90 @@ -421,6 +421,7 @@ module constants ! indicates that an array index hasn't been set integer, parameter :: NONE = 0 + integer, parameter :: C_NONE = -1 ! Codes for read errors -- better hope these numbers are never used in an ! input file! diff --git a/src/geometry.F90 b/src/geometry.F90 index 7fdd03ad8..44f138e10 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -222,7 +222,7 @@ contains ! Move particle to next level and set universe j = j + 1 p % n_coord = j - p % coord(j) % universe = c % fill + p % coord(j) % universe = c % fill() + 1 ! Apply translation if (allocated(c % translation)) then @@ -243,7 +243,7 @@ contains ! ====================================================================== ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL - associate (lat => lattices(c % fill) % obj) + 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) @@ -252,7 +252,7 @@ contains p % coord(j + 1) % uvw = p % coord(j) % uvw ! set particle lattice indices - p % coord(j + 1) % lattice = c % fill + 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) diff --git a/src/geometry_aux.cpp b/src/geometry_aux.cpp index 6efe1cf99..2bf379a6d 100644 --- a/src/geometry_aux.cpp +++ b/src/geometry_aux.cpp @@ -27,10 +27,10 @@ adjust_indices_c() auto search_lat = lattice_dict.find(id); if (search_univ != universe_dict.end()) { c->type = FILL_UNIVERSE; - c->fill = search_univ->second + 1; //TODO: off-by-one + c->fill = search_univ->second; } else if (search_lat != lattice_dict.end()) { c->type = FILL_LATTICE; - c->fill = search_lat->second + 1; //TODO: off-by-one + c->fill = search_lat->second; } else { std::stringstream err_msg; err_msg << "Specified fill " << id << " on cell " << c->id @@ -133,11 +133,11 @@ count_cell_instances(int32_t univ_indx) if (c.type == FILL_UNIVERSE) { // This cell contains another universe. Recurse into that universe. - count_cell_instances(c.fill-1); // TODO: off-by-one + count_cell_instances(c.fill); } else if (c.type == FILL_LATTICE) { // This cell contains a lattice. Recurse into the lattice universes. - Lattice &lat = *lattices_c[c.fill-1]; // TODO: off-by-one + Lattice &lat = *lattices_c[c.fill]; for (auto it = lat.begin(); it != lat.end(); ++it) { count_cell_instances(*it); } @@ -160,11 +160,11 @@ count_universe_instances(int32_t search_univ, int32_t target_univ_id) Cell &c = *cells_c[cell_indx]; if (c.type == FILL_UNIVERSE) { - int32_t next_univ = c.fill - 1; // TODO: off-by-one + int32_t next_univ = c.fill; count += count_universe_instances(next_univ, target_univ_id); } else if (c.type == FILL_LATTICE) { - Lattice &lat = *lattices_c[c.fill - 1]; //TODO: off-by-one + Lattice &lat = *lattices_c[c.fill]; for (auto it = lat.begin(); it != lat.end(); ++it) { int32_t next_univ = *it; count += count_universe_instances(next_univ, target_univ_id); @@ -187,11 +187,11 @@ fill_offset_tables(int32_t target_univ_id, int map) if (c.type == FILL_UNIVERSE) { c.offset[map] = offset; - int32_t search_univ = c.fill - 1; // TODO: off-by-one + int32_t search_univ = c.fill; offset += count_universe_instances(search_univ, target_univ_id); } else if (c.type == FILL_LATTICE) { - Lattice &lat = *lattices_c[c.fill - 1]; // TODO: off-by-one + Lattice &lat = *lattices_c[c.fill]; offset = lat.fill_offset_table(offset, target_univ_id, map); } } @@ -208,10 +208,10 @@ maximum_levels(int32_t univ) for (int32_t cell_indx : universes_c[univ]->cells) { Cell &c = *cells_c[cell_indx]; if (c.type == FILL_UNIVERSE) { - int32_t next_univ = c.fill - 1; // TODO: off-by-one + int32_t next_univ = c.fill; levels_below = std::max(levels_below, maximum_levels(next_univ)); } else if (c.type == FILL_LATTICE) { - Lattice &lat = *lattices_c[c.fill - 1]; //TODO: off-by-one + Lattice &lat = *lattices_c[c.fill]; for (auto it = lat.begin(); it != lat.end(); ++it) { int32_t next_univ = *it; levels_below = std::max(levels_below, maximum_levels(next_univ)); diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index fc70a8eec..d019e8fce 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -18,35 +18,30 @@ module geometry_header interface function cell_pointer_c(cell_ind) bind(C, name='cell_pointer') result(ptr) import C_PTR, C_INT32_T - implicit none integer(C_INT32_T), intent(in), value :: cell_ind type(C_PTR) :: ptr end function cell_pointer_c function cell_id_c(cell_ptr) bind(C, name='cell_id') result(id) import C_PTR, C_INT32_T - implicit none type(C_PTR), intent(in), value :: cell_ptr integer(C_INT32_T) :: id end function cell_id_c subroutine cell_set_id_c(cell_ptr, id) bind(C, name='cell_set_id') import C_PTR, C_INT32_T - implicit none type(C_PTR), intent(in), value :: cell_ptr integer(C_INT32_T), intent(in), value :: id end subroutine cell_set_id_c function cell_type_c(cell_ptr) bind(C, name='cell_type') result(type) import C_PTR, C_INT - implicit none type(C_PTR), intent(in), value :: cell_ptr integer(C_INT) :: type end function cell_type_c subroutine cell_set_type_c(cell_ptr, type) bind(C, name='cell_set_type') import C_PTR, C_INT - implicit none type(C_PTR), intent(in), value :: cell_ptr integer(C_INT), intent(in), value :: type end subroutine cell_set_type_c @@ -54,7 +49,6 @@ module geometry_header function cell_universe_c(cell_ptr) bind(C, name='cell_universe') & result(universe) import C_PTR, C_INT32_T - implicit none type(C_PTR), intent(in), value :: cell_ptr integer(C_INT32_T) :: universe end function cell_universe_c @@ -62,22 +56,31 @@ module geometry_header subroutine cell_set_universe_c(cell_ptr, universe) & bind(C, name='cell_set_universe') import C_PTR, C_INT32_T - implicit none type(C_PTR), intent(in), value :: cell_ptr integer(C_INT32_T), intent(in), value :: universe end subroutine cell_set_universe_c + function cell_fill_c(cell_ptr) bind(C, name="cell_fill") result(fill) + import C_PTR, C_INT32_T + type(C_PTR), intent(in), value :: cell_ptr + integer(C_INT32_T) :: fill + end function cell_fill_c + + function cell_fill_ptr(cell_ptr) bind(C) result(fill_ptr) + import C_PTR + type(C_PTR), intent(in), value :: cell_ptr + type(C_PTR) :: fill_ptr + end function cell_fill_ptr + function cell_n_instances_c(cell_ptr) bind(C, name='cell_n_instances') & result(n_instances) import C_PTR, C_INT32_T - implicit none type(C_PTR), intent(in), value :: cell_ptr integer(C_INT32_T) :: n_instances end function cell_n_instances_c function cell_simple_c(cell_ptr) bind(C, name='cell_simple') result(simple) import C_PTR, C_BOOL - implicit none type(C_PTR), intent(in), value :: cell_ptr logical(C_BOOL) :: simple end function cell_simple_c @@ -85,7 +88,6 @@ module geometry_header subroutine cell_distance_c(cell_ptr, xyz, uvw, on_surface, min_dist, & i_surf) bind(C, name="cell_distance") import C_PTR, C_INT32_T, C_DOUBLE - implicit none type(C_PTR), intent(in), value :: cell_ptr real(C_DOUBLE), intent(in) :: xyz(3) real(C_DOUBLE), intent(in) :: uvw(3) @@ -259,7 +261,6 @@ module geometry_header type Cell type(C_PTR) :: ptr - integer :: fill ! universe # filling this cell integer, allocatable :: material(:) ! Material within cell. Multiple ! materials for distribcell ! instances. 0 signifies a universe @@ -284,6 +285,7 @@ module geometry_header procedure :: set_type => cell_set_type procedure :: universe => cell_universe procedure :: set_universe => cell_set_universe + procedure :: fill => cell_fill procedure :: n_instances => cell_n_instances procedure :: simple => cell_simple procedure :: distance => cell_distance @@ -414,6 +416,12 @@ contains call cell_set_universe_c(this % ptr, universe) end subroutine cell_set_universe + function cell_fill(this) result(fill) + class(Cell), intent(in) :: this + integer(C_INT32_T) :: fill + fill = cell_fill_c(this % ptr) + end function cell_fill + function cell_n_instances(this) result(n_instances) class(Cell), intent(in) :: this integer(C_INT32_T) :: n_instances @@ -611,7 +619,7 @@ contains indices = C_LOC(c % material(1)) case (FILL_UNIVERSE, FILL_LATTICE) n = 1 - indices = C_LOC(c % fill) + indices = cell_fill_ptr(c % ptr) end select end associate else @@ -723,7 +731,7 @@ contains if (index >= 1 .and. index <= size(cells)) then ! error if the cell is filled with another universe - if (cells(index) % fill /= NONE) then + if (cells(index) % fill() /= C_NONE) then err = E_GEOMETRY call set_errmsg("Cannot set temperature on a cell filled & &with a universe.") diff --git a/src/input_xml.F90 b/src/input_xml.F90 index d325928a6..501a43a57 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -1056,12 +1056,6 @@ contains ! Get pointer to i-th cell node node_cell = node_cell_list(i) - if (check_for_node(node_cell, "fill")) then - call get_node_value(node_cell, "fill", c % fill) - else - c % fill = NONE - end if - ! Check to make sure 'id' hasn't been used if (cell_dict % has(c % id())) then call fatal_error("Two or more cells use the same unique ID: " & @@ -1104,18 +1098,6 @@ contains c % material(1) = NONE end if - ! Check to make sure that either material or fill was specified - if (c % material(1) == NONE .and. c % fill == NONE) then - call fatal_error("Neither material nor fill was specified for cell " & - // trim(to_str(c % id()))) - end if - - ! Check to make sure that both material and fill haven't been - ! specified simultaneously - if (c % material(1) /= NONE .and. c % fill /= NONE) then - call fatal_error("Cannot specify material and fill simultaneously") - end if - ! Check for region specification (also under deprecated name surfaces) if (check_for_node(node_cell, "surfaces")) then call warning("The use of 'surfaces' is deprecated and will be & @@ -1157,7 +1139,7 @@ contains if (check_for_node(node_cell, "rotation")) then ! Rotations can only be applied to cells that are being filled with ! another universe - if (c % fill == NONE) then + if (c % fill() == C_NONE) then call fatal_error("Cannot apply a rotation to cell " // trim(to_str(& &c % id())) // " because it is not filled with another universe") end if @@ -1192,7 +1174,7 @@ contains if (check_for_node(node_cell, "translation")) then ! Translations can only be applied to cells that are being filled with ! another universe - if (c % fill == NONE) then + if (c % fill() == C_NONE) then call fatal_error("Cannot apply a translation to cell " & // trim(to_str(c % id())) // " because it is not filled with & &another universe") @@ -3800,7 +3782,6 @@ contains integer :: i ! index for various purposes integer :: j ! index for various purposes - integer :: lid ! lattice IDs integer :: id ! user-specified id call adjust_indices_c() @@ -3811,15 +3792,7 @@ contains ! ======================================================================= ! ADJUST MATERIAL/FILL POINTERS FOR EACH CELL - if (c % material(1) == NONE) then - id = c % fill - if (universe_dict % has(id)) then - c % fill = universe_dict % get(id) - elseif (lattice_dict % has(id)) then - lid = lattice_dict % get(id) - c % fill = lid - end if - else + if (c % material(1) /= NONE) then do j = 1, size(c % material) id = c % material(j) if (id == MATERIAL_VOID) then diff --git a/src/summary.F90 b/src/summary.F90 index 1d3349123..3b01ddfc9 100644 --- a/src/summary.F90 +++ b/src/summary.F90 @@ -160,7 +160,7 @@ contains subroutine write_geometry(file_id) integer(HID_T), intent(in) :: file_id - integer :: i, j + integer :: i, j, cell_fill integer, allocatable :: cell_materials(:) integer, allocatable :: cell_ids(:) real(8), allocatable :: cell_temperatures(:) @@ -225,7 +225,7 @@ contains case (FILL_UNIVERSE) call write_dataset(cell_group, "fill_type", "universe") - call write_dataset(cell_group, "fill", universes(c%fill)%id) + call write_dataset(cell_group, "fill", universes(c%fill()+1)%id) if (allocated(c%translation)) then call write_dataset(cell_group, "translation", c%translation) @@ -236,7 +236,10 @@ contains case (FILL_LATTICE) call write_dataset(cell_group, "fill_type", "lattice") - call write_dataset(cell_group, "lattice", lattices(c%fill)%obj%id()) + ! Do not access the 'lattices' array with 'c % fill() + 1' directly; it + ! causes a segfault in GCC 7.3.0 + cell_fill = c % fill() + 1 + call write_dataset(cell_group, "lattice", lattices(cell_fill)%obj%id()) end select call close_group(cell_group) diff --git a/src/tallies/tally_filter_distribcell.F90 b/src/tallies/tally_filter_distribcell.F90 index ce2549f58..fa0a6b296 100644 --- a/src/tallies/tally_filter_distribcell.F90 +++ b/src/tallies/tally_filter_distribcell.F90 @@ -215,7 +215,7 @@ contains ! Get the offset of the first lattice location else - lat => lattices(c % fill) % obj + lat => lattices(c % fill() + 1) % obj temp_offset = lat % offset(map-1, [0, 0, 0]) end if @@ -250,7 +250,7 @@ contains ! Enter this cell to update the current offset offset = c % offset(map-1) + offset - next_univ => universes(c % fill) + next_univ => universes(c % fill() + 1) call find_offset(i_cell, next_univ, target_offset, offset, path) return @@ -259,7 +259,7 @@ contains elseif (c % type() == FILL_LATTICE) then ! Set current lattice - lat => lattices(c % fill) % obj + lat => lattices(c % fill() + 1) % obj select type (lat)