diff --git a/src/geometry.F90 b/src/geometry.F90 index ac9f74bb34..54557dfd16 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -2,7 +2,8 @@ module geometry use constants use error, only: fatal_error, warning - use geometry_header, only: Cell, Surface, Universe, Lattice + use geometry_header, only: Cell, Surface, Universe, Lattice, & + &RectLattice, HexLattice use global use output, only: write_message use particle_header, only: LocalCoord, deallocate_coord, Particle @@ -138,7 +139,7 @@ contains real(8) :: xyz(3) ! temporary location logical :: use_search_cells ! use cells provided as argument type(Cell), pointer, save :: c => null() ! pointer to cell - type(Lattice), pointer, save :: lat => null() ! pointer to lattice + class(Lattice), pointer, save :: lat => null() ! pointer to lattice type(Universe), pointer, save :: univ => null() ! universe to search in !$omp threadprivate(c, lat, univ) @@ -222,7 +223,7 @@ contains ! CELL CONTAINS LATTICE, RECURSIVELY FIND CELL ! Set current lattice - lat => lattices(c % fill) + lat => lattices(c % fill) % obj ! Determine lattice indices i_xyz = get_lat_indices(lat, & @@ -279,15 +280,17 @@ contains !=============================================================================== function get_lat_trans(lat, xyz, i_xyz) result(xyz_t) - type(Lattice), pointer, intent(in) :: lat - real(8) , intent(in) :: xyz(3) - integer , intent(in) :: i_xyz(3) + class(Lattice), pointer, intent(in) :: lat + real(8) , intent(in) :: xyz(3) + integer , intent(in) :: i_xyz(3) real(8) :: xyz_t(3) integer :: n_rings - if (lat % type == LATTICE_RECT) then + select type(lat) + + type is (RectLattice) xyz_t(1) = xyz(1) - (lat % lower_left(1) + & &(i_xyz(1) - 0.5_8)*lat % width(1)) xyz_t(2) = xyz(2) - (lat % lower_left(2) + & @@ -299,7 +302,7 @@ contains xyz_t(3) = xyz(3) end if - else if (lat % type == LATTICE_HEX) then + type is (HexLattice) n_rings = lat % dimension(1) xyz_t(1) = xyz(1) - (lat % lower_left(1) + & @@ -314,7 +317,7 @@ contains xyz_t(3) = xyz(3) end if - end if + end select end function get_lat_trans @@ -324,13 +327,15 @@ contains !=============================================================================== function is_valid_lat_index(lat, i_xyz) result(is_valid) - type(Lattice), pointer, intent(in) :: lat - integer , intent(in) :: i_xyz(3) - logical :: is_valid + class(Lattice), pointer, intent(in) :: lat + integer , intent(in) :: i_xyz(3) + logical :: is_valid integer :: n_rings, n_x, n_y, n_z - if (lat % type == LATTICE_RECT) then + select type(lat) + + type is (RectLattice) n_x = lat % dimension(1) n_y = lat % dimension(2) if (lat % n_dimension == 3) then @@ -343,7 +348,7 @@ contains &i_xyz(2) > 0 .and. i_xyz(2) <= n_y .and. & &i_xyz(3) > 0 .and. i_xyz(3) <= n_z - else if (lat % type == LATTICE_HEX) then + type is (HexLattice) n_rings = lat % dimension(1) if (lat % n_dimension == 2) then @@ -358,7 +363,7 @@ contains &i_xyz(1) + i_xyz(2) < 3*n_rings .and. & &i_xyz(3) > 0 .and. i_xyz(3) < n_z + 1) - end if + end select end function is_valid_lat_index @@ -368,9 +373,9 @@ contains !=============================================================================== function get_lat_indices(lat, xyz, uvw) result(i_xyz) - type(Lattice), pointer, intent(in) :: lat - real(8) , intent(in) :: xyz(3) - real(8) , intent(in) :: uvw(3) + class(Lattice), pointer, intent(in) :: lat + real(8) , intent(in) :: xyz(3) + real(8) , intent(in) :: uvw(3) integer :: i_xyz(3) @@ -379,8 +384,10 @@ contains real(8) :: dists(4) integer :: n_rings integer :: i, j, k, loc(1) - - if (lat % type == LATTICE_RECT) then + + select type(lat) + + type is (RectLattice) ! Find approximate indices using ceiling division. i_xyz(1) = ceiling((xyz(1) - lat % lower_left(1))/lat % width(1)) i_xyz(2) = ceiling((xyz(2) - lat % lower_left(2))/lat % width(2)) @@ -416,7 +423,7 @@ contains ! end if !end if - else if (lat % type == LATTICE_HEX) then + type is (HexLattice) ! Index z direction. if (lat % n_dimension == 2) then i_xyz(3) = ceiling((xyz(3) - lat % lower_left(3)) / lat % width(2)) @@ -533,7 +540,7 @@ contains end if - end if + end select end function get_lat_indices @@ -849,11 +856,11 @@ contains real(8) :: x0, y0, z0 ! half width of lattice element real(8) :: half_pitch ! half pitch of a hexagonal lattice logical :: found ! particle found in cell? - type(Lattice), pointer, save :: lat => null() + class(Lattice), pointer, save :: lat => null() type(LocalCoord), pointer :: parent_coord !$omp threadprivate(lat) - lat => lattices(p % coord % lattice) + lat => lattices(p % coord % lattice) % obj if (verbosity >= 10 .or. trace) then message = " Crossing lattice " // trim(to_str(lat % id)) // & @@ -952,7 +959,7 @@ contains logical :: on_surface ! is particle on surface? type(Cell), pointer, save :: cl => null() type(Surface), pointer, save :: surf => null() - type(Lattice), pointer, save :: lat => null() + class(Lattice), pointer, save :: lat => null() type(LocalCoord), pointer, save :: coord => null() type(LocalCoord), pointer, save :: final_coord => null() @@ -1413,8 +1420,11 @@ contains ! FIND MINIMUM DISTANCE TO LATTICE SURFACES LAT_COORD: if (coord % lattice /= NONE) then - lat => lattices(coord % lattice) - LAT_TYPE: if (lat % type == LATTICE_RECT) then + lat => lattices(coord % lattice) % obj + + LAT_TYPE: select type(lat) + + type is (RectLattice) ! copy local coordinates x = coord % xyz(1) y = coord % xyz(2) @@ -1480,7 +1490,7 @@ contains end if end if - elseif (lat % type == LATTICE_HEX) then LAT_TYPE + type is (HexLattice) LAT_TYPE ! Copy local coordinates. x = coord % xyz(1) y = coord % xyz(2) @@ -1592,7 +1602,7 @@ contains end if end if end if - end if LAT_TYPE + end select LAT_TYPE ! If the lattice boundary is coincident with the parent cell boundary, ! we need to make sure that the lattice is not selected. This is diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 93863865c4..7ac32cae73 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -21,18 +21,28 @@ module geometry_header ! triangular) !=============================================================================== - type Lattice - integer :: id ! Universe number for lattice - integer :: type ! Type of lattice (rectangular, hex, etc) - integer :: level ! Level of lattice - integer :: n_dimension ! Number of dimensions - integer, allocatable :: dimension(:) ! number of cells in each direction - real(8), allocatable :: lower_left(:) ! lower-left corner of lattice - real(8), allocatable :: width(:) ! width of each lattice cell - integer, allocatable :: universes(:,:,:) ! specified universes - integer :: outside ! material to fill area outside + type, abstract :: Lattice + integer :: id ! Universe number for lattice + integer :: type ! Type of lattice (rectangular, hex, etc) + integer :: level ! Level of lattice + integer :: n_dimension ! Number of dimensions + integer, allocatable :: dimension(:) ! number of cells in each direction + real(8), allocatable :: lower_left(:) ! lower-left corner of lattice + real(8), allocatable :: width(:) ! width of each lattice cell + integer, allocatable :: universes(:,:,:) ! specified universes + integer :: outside ! material to fill area outside end type Lattice + type, extends(Lattice) :: RectLattice + end type RectLattice + + type, extends(Lattice) :: HexLattice + end type HexLattice + + type LatticeContainer + class(Lattice), allocatable :: obj + end type LatticeContainer + !=============================================================================== ! SURFACE type defines a first- or second-order surface that can be used to ! construct closed volumes (cells) diff --git a/src/global.F90 b/src/global.F90 index 7d728a25d6..4a21ffe4d4 100644 --- a/src/global.F90 +++ b/src/global.F90 @@ -6,7 +6,7 @@ module global use cmfd_header use constants use dict_header, only: DictCharInt, DictIntInt - use geometry_header, only: Cell, Universe, Lattice, Surface + use geometry_header, only: Cell, Universe, Lattice, LatticeContainer, Surface use material_header, only: Material use mesh_header, only: StructuredMesh use plot_header, only: ObjectPlot @@ -26,12 +26,12 @@ module global ! GEOMETRY-RELATED VARIABLES ! Main arrays - type(Cell), allocatable, target :: cells(:) - type(Universe), allocatable, target :: universes(:) - type(Lattice), allocatable, target :: lattices(:) - type(Surface), allocatable, target :: surfaces(:) - type(Material), allocatable, target :: materials(:) - type(ObjectPlot),allocatable, target :: plots(:) + type(Cell), allocatable, target :: cells(:) + type(Universe), allocatable, target :: universes(:) + type(LatticeContainer), allocatable, target :: lattices(:) + type(Surface), allocatable, target :: surfaces(:) + type(Material), allocatable, target :: materials(:) + type(ObjectPlot), allocatable, target :: plots(:) ! Size of main arrays integer :: n_cells ! # of cells diff --git a/src/hdf5_summary.F90 b/src/hdf5_summary.F90 index c9e367200c..adaf5b2412 100644 --- a/src/hdf5_summary.F90 +++ b/src/hdf5_summary.F90 @@ -105,7 +105,7 @@ contains type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() type(Universe), pointer :: u => null() - type(Lattice), pointer :: lat => null() + class(Lattice), pointer :: lat => null() ! Use H5LT interface to write number of geometry objects call su % write_data(n_cells, "n_cells", group="geometry") @@ -275,11 +275,11 @@ contains lat => lattices(i) ! Write lattice type - select case(lat % type) - case (LATTICE_RECT) + select type (lat) + type is (LatticeRect) call su % write_data("rectangular", "type", & group="geometry/lattices/lattice " // trim(to_str(lat % id))) - case (LATTICE_HEX) + type is (LatticeHex) call su % write_data("hexagonal", "type", & group="geometry/lattices/lattice " // trim(to_str(lat % id))) end select diff --git a/src/initialize.F90 b/src/initialize.F90 index acb654bb98..6ab3ac4930 100644 --- a/src/initialize.F90 +++ b/src/initialize.F90 @@ -7,7 +7,8 @@ module initialize use energy_grid, only: unionized_grid use error, only: fatal_error, warning use geometry, only: neighbor_lists - use geometry_header, only: Cell, Universe, Lattice, BASE_UNIVERSE + use geometry_header, only: Cell, Universe, Lattice, RectLattice, HexLattice,& + &BASE_UNIVERSE use global use input_xml, only: read_input_xml, read_cross_sections_xml, & cells_in_univ_dict, read_plots_xml @@ -566,7 +567,7 @@ contains integer :: i_array ! index in surfaces/materials array integer :: id ! user-specified id type(Cell), pointer :: c => null() - type(Lattice), pointer :: lat => null() + class(Lattice), pointer :: lat => null() type(TallyObject), pointer :: t => null() do i = 1, n_cells @@ -622,7 +623,7 @@ contains c % fill = universe_dict % get_key(id) elseif (lattice_dict % has_key(id)) then lid = lattice_dict % get_key(id) - mid = lattices(lid) % outside + mid = lattices(lid) % obj % outside c % type = CELL_LATTICE c % fill = lid if (mid == MATERIAL_VOID) then @@ -647,9 +648,10 @@ contains ! ADJUST UNIVERSE INDICES FOR EACH LATTICE do i = 1, n_lattices - lat => lattices(i) + lat => lattices(i) % obj + select type (lat) - if (lat % type == LATTICE_RECT) then + type is (RectLattice) n_x = lat % dimension(1) n_y = lat % dimension(2) if (lat % n_dimension == 3) then @@ -673,7 +675,7 @@ contains end do end do - else + type is (HexLattice) n_rings = lat % dimension(1) if (lat % n_dimension == 2) then n_z = lat % dimension(2) @@ -701,7 +703,7 @@ contains end do end do - end if + end select end do diff --git a/src/input_xml.F90 b/src/input_xml.F90 index 4924346549..366e5f1a72 100644 --- a/src/input_xml.F90 +++ b/src/input_xml.F90 @@ -4,7 +4,7 @@ module input_xml use constants use dict_header, only: DictIntInt, ElemKeyValueCI use error, only: fatal_error, warning - use geometry_header, only: Cell, Surface, Lattice + use geometry_header, only: Cell, Surface, Lattice, RectLattice, HexLattice use global use list_header, only: ListChar, ListReal use mesh_header, only: StructuredMesh @@ -827,7 +827,7 @@ contains integer :: i, j, k, m, i_x, i_a, input_index integer :: n - integer :: n_x, n_y, n_z, n_rings + integer :: n_x, n_y, n_z, n_rings, n_rlats, n_hlats integer :: universe_num integer :: n_cells_in_univ integer :: coeffs_reqd @@ -839,16 +839,17 @@ contains logical :: boundary_exists character(MAX_LINE_LEN) :: filename character(MAX_WORD_LEN) :: word - type(Cell), pointer :: c => null() - type(Surface), pointer :: s => null() - type(Lattice), pointer :: lat => null() + type(Cell), pointer :: c => null() + type(Surface), pointer :: s => null() + class(Lattice), pointer :: lat => null() type(Node), pointer :: doc => null() type(Node), pointer :: node_cell => null() type(Node), pointer :: node_surf => null() type(Node), pointer :: node_lat => null() type(NodeList), pointer :: node_cell_list => null() type(NodeList), pointer :: node_surf_list => null() - type(NodeList), pointer :: node_lat_list => null() + type(NodeList), pointer :: node_rlat_list => null() + type(NodeList), pointer :: node_hlat_list => null() ! Display output message message = "Reading geometry XML file..." @@ -1187,17 +1188,21 @@ contains ! READ LATTICES FROM GEOMETRY.XML ! Get pointer to list of XML - call get_node_list(doc, "lattice", node_lat_list) + call get_node_list(doc, "lattice", node_rlat_list) + call get_node_list(doc, "hex_lattice", node_hlat_list) ! Allocate lattices array - n_lattices = get_list_size(node_lat_list) + n_rlats = get_list_size(node_rlat_list) + n_hlats = get_list_size(node_hlat_list) + n_lattices = n_rlats + n_hlats allocate(lattices(n_lattices)) - do i = 1, n_lattices - lat => lattices(i) + RECT_LATTICES: do i = 1, n_rlats + allocate(RectLattice::lattices(i) % obj) + lat => lattices(i) % obj ! Get pointer to i-th lattice - call get_list_item(node_lat_list, i, node_lat) + call get_list_item(node_rlat_list, i, node_lat) ! ID of lattice if (check_for_node(node_lat, "id")) then @@ -1214,28 +1219,11 @@ contains call fatal_error() end if - ! Read lattice type - word = '' - if (check_for_node(node_lat, "type")) & - call get_node_value(node_lat, "type", word) - call lower_case(word) - select case (trim(word)) - case ('rect', 'rectangle', 'rectangular') - lat % type = LATTICE_RECT - case ('hex', 'hexagon', 'hexagonal') - lat % type = LATTICE_HEX - case default - message = "Invalid lattice type: " // trim(word) - call fatal_error() - end select - ! Read number of lattice cells in each dimension n = get_arraysize_integer(node_lat, "dimension") - if (lat % type == LATTICE_RECT .and. n /= 2 .and. n /= 3) then + if (n /= 2 .and. n /= 3) then message = "Rectangular lattice must be two or three dimensions." call fatal_error() - else if (lat % type == LATTICE_HEX .and. n /= 1 .and. n/= 2) then - message = "Hexagonal lattice must be one or two dimension(s)." end if lat % n_dimension = n @@ -1243,28 +1231,16 @@ contains call get_node_array(node_lat, "dimension", lat % dimension) ! Read lattice lower-left location - if (lat % type == LATTICE_RECT) then - if (size(lat % dimension) /= & - &get_arraysize_double(node_lat, "lower_left")) then - message = "Number of entries on must be the same as & - &the number of entries on ." - call fatal_error() - end if - - allocate(lat % lower_left(n)) - call get_node_array(node_lat, "lower_left", lat % lower_left) - else - if (size(lat % dimension) /= & - &get_arraysize_double(node_lat, "lower_left") - 1) then - message = "Number of entries on must be one greater & - &than the number of entries on ." - call fatal_error() - end if - - allocate(lat % lower_left(n+1)) - call get_node_array(node_lat, "lower_left", lat % lower_left) + if (size(lat % dimension) /= & + &get_arraysize_double(node_lat, "lower_left")) then + message = "Number of entries on must be the same as & + &the number of entries on ." + call fatal_error() end if + allocate(lat % lower_left(n)) + call get_node_array(node_lat, "lower_left", lat % lower_left) + ! Read lattice widths if (size(lat % dimension) /= & get_arraysize_double(node_lat, "width")) then @@ -1277,145 +1253,35 @@ contains call get_node_array(node_lat, "width", lat % width) ! Copy number of dimensions - if (lat % type == LATTICE_RECT) then - n_x = lat % dimension(1) - n_y = lat % dimension(2) - if (lat % n_dimension == 3) then - n_z = lat % dimension(3) - else - n_z = 1 - end if - allocate(lat % universes(n_x, n_y, n_z)) + n_x = lat % dimension(1) + n_y = lat % dimension(2) + if (lat % n_dimension == 3) then + n_z = lat % dimension(3) else - n_rings = lat % dimension(1) - if (lat % n_dimension == 2) then - n_z = lat % dimension(2) - else - n_z = 1 - end if - allocate(lat % universes(2*n_rings - 1, 2*n_rings - 1, n_z)) + n_z = 1 end if + allocate(lat % universes(n_x, n_y, n_z)) ! Check that number of universes matches size n = get_arraysize_integer(node_lat, "universes") - if (lat % type == LATTICE_RECT) then - if (n /= n_x*n_y*n_z) then - message = "Number of universes on does not match size of & - &lattice " // trim(to_str(lat % id)) // "." - call fatal_error() - end if - else - if (n /= (3*n_rings**2 - 3*n_rings + 1)*n_z) then - message = "Number of universes on does not match size of & - &lattice " // trim(to_str(lat % id)) // "." - call fatal_error() - end if + if (n /= n_x*n_y*n_z) then + message = "Number of universes on does not match size of & + &lattice " // trim(to_str(lat % id)) // "." + call fatal_error() end if allocate(temp_int_array(n)) call get_node_array(node_lat, "universes", temp_int_array) ! Read universes - if (lat % type == LATTICE_RECT) then - do m = 1, n_z - do k = 0, n_y - 1 - do j = 1, n_x - lat % universes(j, n_y - k, m) = & - temp_int_array(j + n_x*k + n_x*n_y*(m-1)) - end do + do m = 1, n_z + do k = 0, n_y - 1 + do j = 1, n_x + lat % universes(j, n_y - k, m) = & + &temp_int_array(j + n_x*k + n_x*n_y*(m-1)) end do end do - else - do m = 1, n_z - do k = 1, 2*n_rings - 1 - do j = 1, 2*n_rings - 1 - lat % universes(j, k, m) = -1 - end do - end do - end do - ! TODO: The index walk in the k loops can be made a little more - ! efficient. They currently sometimes change index values and then - ! change them again before use. - - ! Universes in hexagonal lattices are stored in a manner that represents - ! a skewed coordinate system: (x, alpha) rather than (x, y). There is - ! no obvious, direct relationship between the order of universes in the - ! input and the order that they will be stored in the skewed array so - ! the following code walks a set of index values across the skewed array - ! in a manner that matches the input order. Note that i_x = 0, i_a = 0 - ! corresponds to the center of the hexagonal lattice. - - input_index = 1 - do m = 1, n_z - ! Initialize lattice indecies. - i_x = 1 - i_a = n_rings - 1 - - ! Map upper triangular region of hexagonal lattice. - do k = 1, n_rings-1 - ! Walk index to lower-left neighbor of last row start. - i_x = i_x - 1 - do j = 1, k - ! Place universe in array. - lat % universes(i_x + n_rings, i_a + n_rings, m) = & - &temp_int_array(input_index) - ! Walk index to closest non-adjacent right neighbor. - i_x = i_x + 2 - i_a = i_a - 1 - ! Increment XML array index. - input_index = input_index + 1 - end do - ! Return lattice index to start of current row. - i_x = i_x - 2*k - i_a = i_a + k - end do - - ! Map middle square region of hexagonal lattice. - do k = 1, 2*n_rings - 1 - if (mod(k, 2) == 1) then - ! Walk index to lower-left neighbor of last row start. - i_x = i_x - 1 - else - ! Walk index to lower-right neighbor of last row start - i_x = i_x + 1 - i_a = i_a - 1 - end if - do j = 1, n_rings - mod(k-1, 2) - ! Place universe in array. - lat % universes(i_x + n_rings, i_a + n_rings, m) = & - &temp_int_array(input_index) - ! Walk index to closest non-adjacent right neighbor. - i_x = i_x + 2 - i_a = i_a - 1 - ! Increment XML array index. - input_index = input_index + 1 - end do - ! Return lattice index to start of current row. - i_x = i_x - 2*(n_rings - mod(k-1, 2)) - i_a = i_a + n_rings - mod(k-1, 2) - end do - - ! Map lower triangular region of hexagonal lattice. - do k = 1, n_rings-1 - ! Walk index to lower-right neighbor of last row start. - i_x = i_x + 1 - i_a = i_a - 1 - do j = 1, n_rings - k - ! Place universe in array. - lat % universes(i_x + n_rings, i_a + n_rings, m) = & - &temp_int_array(input_index) - ! Walk index to closest non-adjacent right neighbor. - i_x = i_x + 2 - i_a = i_a - 1 - ! Increment XML array index. - input_index = input_index + 1 - end do - ! Return lattice index to start of current row. - i_x = i_x - 2*(n_rings - k) - i_a = i_a + n_rings - k - end do - end do - end if + end do deallocate(temp_int_array) ! Read material for area outside lattice @@ -1432,7 +1298,183 @@ contains ! Add lattice to dictionary call lattice_dict % add_key(lat % id, i) - end do + end do RECT_LATTICES + + HEX_LATTICES: do i = 1, n_hlats + allocate(HexLattice::lattices(i) % obj) + lat => lattices(i) % obj + + ! Get pointer to i-th lattice + call get_list_item(node_hlat_list, i, node_lat) + + ! ID of lattice + if (check_for_node(node_lat, "id")) then + call get_node_value(node_lat, "id", lat % id) + else + message = "Must specify id of lattice in geometry XML file." + call fatal_error() + end if + + ! Check to make sure 'id' hasn't been used + if (lattice_dict % has_key(lat % id)) then + message = "Two or more lattices use the same unique ID: " // & + to_str(lat % id) + call fatal_error() + end if + + ! Read number of lattice cells in each dimension + n = get_arraysize_integer(node_lat, "dimension") + if (n /= 1 .and. n/= 2) then + message = "Hexagonal lattice must be one or two dimension(s)." + call fatal_error() + end if + + lat % n_dimension = n + allocate(lat % dimension(n)) + call get_node_array(node_lat, "dimension", lat % dimension) + + ! Read lattice lower-left location + if (size(lat % dimension) /= & + &get_arraysize_double(node_lat, "lower_left") - 1) then + message = "Number of entries on must be one greater & + &than the number of entries on ." + call fatal_error() + end if + + allocate(lat % lower_left(n+1)) + call get_node_array(node_lat, "lower_left", lat % lower_left) + + ! Read lattice widths + if (size(lat % dimension) /= & + get_arraysize_double(node_lat, "width")) then + message = "Number of entries on must be the same as & + &the number of entries on ." + call fatal_error() + end if + + allocate(lat % width(n)) + call get_node_array(node_lat, "width", lat % width) + + ! Copy number of dimensions + n_rings = lat % dimension(1) + if (lat % n_dimension == 2) then + n_z = lat % dimension(2) + else + n_z = 1 + end if + allocate(lat % universes(2*n_rings - 1, 2*n_rings - 1, n_z)) + + ! Check that number of universes matches size + n = get_arraysize_integer(node_lat, "universes") + if (n /= (3*n_rings**2 - 3*n_rings + 1)*n_z) then + message = "Number of universes on does not match size of & + &lattice " // trim(to_str(lat % id)) // "." + call fatal_error() + end if + + allocate(temp_int_array(n)) + call get_node_array(node_lat, "universes", temp_int_array) + + ! Read universes + ! TODO: The index walk in the k loops can be made a little more + ! efficient. They currently sometimes change index values and then + ! change them again before use. + + ! Universes in hexagonal lattices are stored in a manner that represents + ! a skewed coordinate system: (x, alpha) rather than (x, y). There is + ! no obvious, direct relationship between the order of universes in the + ! input and the order that they will be stored in the skewed array so + ! the following code walks a set of index values across the skewed array + ! in a manner that matches the input order. Note that i_x = 0, i_a = 0 + ! corresponds to the center of the hexagonal lattice. + + input_index = 1 + do m = 1, n_z + ! Initialize lattice indecies. + i_x = 1 + i_a = n_rings - 1 + + ! Map upper triangular region of hexagonal lattice. + do k = 1, n_rings-1 + ! Walk index to lower-left neighbor of last row start. + i_x = i_x - 1 + do j = 1, k + ! Place universe in array. + lat % universes(i_x + n_rings, i_a + n_rings, m) = & + &temp_int_array(input_index) + ! Walk index to closest non-adjacent right neighbor. + i_x = i_x + 2 + i_a = i_a - 1 + ! Increment XML array index. + input_index = input_index + 1 + end do + ! Return lattice index to start of current row. + i_x = i_x - 2*k + i_a = i_a + k + end do + + ! Map middle square region of hexagonal lattice. + do k = 1, 2*n_rings - 1 + if (mod(k, 2) == 1) then + ! Walk index to lower-left neighbor of last row start. + i_x = i_x - 1 + else + ! Walk index to lower-right neighbor of last row start + i_x = i_x + 1 + i_a = i_a - 1 + end if + do j = 1, n_rings - mod(k-1, 2) + ! Place universe in array. + lat % universes(i_x + n_rings, i_a + n_rings, m) = & + &temp_int_array(input_index) + ! Walk index to closest non-adjacent right neighbor. + i_x = i_x + 2 + i_a = i_a - 1 + ! Increment XML array index. + input_index = input_index + 1 + end do + ! Return lattice index to start of current row. + i_x = i_x - 2*(n_rings - mod(k-1, 2)) + i_a = i_a + n_rings - mod(k-1, 2) + end do + + ! Map lower triangular region of hexagonal lattice. + do k = 1, n_rings-1 + ! Walk index to lower-right neighbor of last row start. + i_x = i_x + 1 + i_a = i_a - 1 + do j = 1, n_rings - k + ! Place universe in array. + lat % universes(i_x + n_rings, i_a + n_rings, m) = & + &temp_int_array(input_index) + ! Walk index to closest non-adjacent right neighbor. + i_x = i_x + 2 + i_a = i_a - 1 + ! Increment XML array index. + input_index = input_index + 1 + end do + ! Return lattice index to start of current row. + i_x = i_x - 2*(n_rings - k) + i_a = i_a + n_rings - k + end do + end do + deallocate(temp_int_array) + + ! Read material for area outside lattice + lat % outside = MATERIAL_VOID + if (check_for_node(node_lat, "outside")) then + call get_node_value(node_lat, "outside", mid) + if (mid == 0 .or. mid == MATERIAL_VOID) then + lat % outside = MATERIAL_VOID + else + lat % outside = mid + end if + end if + + ! Add lattice to dictionary + call lattice_dict % add_key(lat % id, i) + + end do HEX_LATTICES ! Close geometry XML file call close_xmldoc(doc) diff --git a/src/output.F90 b/src/output.F90 index 834c28e188..7630228233 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -256,7 +256,7 @@ contains type(Cell), pointer :: c => null() type(Surface), pointer :: s => null() type(Universe), pointer :: u => null() - type(Lattice), pointer :: l => null() + class(Lattice), pointer :: l => null() type(LocalCoord), pointer :: coord => null() ! display type of particle @@ -292,7 +292,7 @@ contains ! Print information on lattice if (coord % lattice /= NONE) then - l => lattices(coord % lattice) + l => lattices(coord % lattice) % obj write(ou,*) ' Lattice = ' // trim(to_str(l % id)) write(ou,*) ' Lattice position = (' // trim(to_str(& p % coord % lattice_x)) // ',' // trim(to_str(& @@ -355,7 +355,7 @@ contains integer :: unit_ ! unit to write to character(MAX_LINE_LEN) :: string type(Universe), pointer :: u => null() - type(Lattice), pointer :: l => null() + class(Lattice), pointer :: l => null() type(Material), pointer :: m => null() ! Set unit to stdout if not already set @@ -384,7 +384,7 @@ contains u => universes(c % fill) write(unit_,*) ' Fill = Universe ' // to_str(u % id) case (CELL_LATTICE) - l => lattices(c % fill) + l => lattices(c % fill) % obj write(unit_,*) ' Fill = Lattice ' // to_str(l % id) end select @@ -471,8 +471,8 @@ contains subroutine print_lattice(lat, unit) - type(Lattice), pointer :: lat - integer, optional :: unit + class(Lattice), pointer :: lat + integer, optional :: unit integer :: i ! loop index integer :: unit_ ! unit to write to @@ -919,7 +919,7 @@ contains type(Surface), pointer :: s => null() type(Cell), pointer :: c => null() type(Universe), pointer :: u => null() - type(Lattice), pointer :: l => null() + class(Lattice), pointer :: l => null() ! print summary of surfaces call header("SURFACE SUMMARY", unit=UNIT_SUMMARY) @@ -946,7 +946,7 @@ contains if (n_lattices > 0) then call header("LATTICE SUMMARY", unit=UNIT_SUMMARY) do i = 1, n_lattices - l => lattices(i) + l => lattices(i) % obj call print_lattice(l, unit=UNIT_SUMMARY) end do end if