From d2d70a3d42764291a68caf7a4c968e2fbe78b6e7 Mon Sep 17 00:00:00 2001 From: Paul Romano Date: Mon, 7 Aug 2017 13:40:49 -0500 Subject: [PATCH] Move most of mesh module into mesh_header as type-bound procedures --- src/cmfd_data.F90 | 7 +- src/cmfd_execute.F90 | 4 +- src/mesh.F90 | 334 +--------------------------------------- src/mesh_header.F90 | 359 +++++++++++++++++++++++++++++++++++++++++++ src/output.F90 | 3 +- src/physics.F90 | 3 +- src/physics_mg.F90 | 3 +- src/plot.F90 | 6 +- src/tally.F90 | 24 +-- src/tally_filter.F90 | 28 +--- src/trigger.F90 | 3 +- 11 files changed, 384 insertions(+), 390 deletions(-) diff --git a/src/cmfd_data.F90 b/src/cmfd_data.F90 index 177f140b63..c873b1beb4 100644 --- a/src/cmfd_data.F90 +++ b/src/cmfd_data.F90 @@ -58,7 +58,6 @@ contains use error, only: fatal_error use global, only: cmfd, n_cmfd_tallies, cmfd_tallies, meshes, & filters, filter_matches - use mesh, only: mesh_indices_to_bin use mesh_header, only: RegularMesh use string, only: to_str use tally_header, only: TallyObject @@ -166,7 +165,7 @@ contains ! Get bin number for mesh indices filter_matches(i_filter_mesh) % bins % data(1) = & - mesh_indices_to_bin(m,ijk) + m % get_bin_from_indices(ijk) ! Apply energy in filter if (energy_filters) then @@ -217,7 +216,7 @@ contains ! Get bin number for mesh indices filter_matches(i_filter_mesh) % bins % data(1) = & - mesh_indices_to_bin(m,ijk) + m % get_bin_from_indices(ijk) if (energy_filters) then ! Apply energy in filter @@ -266,7 +265,7 @@ contains ! Get the bin for this mesh cell filter_matches(i_filter_mesh) % bins % data(1) = & - mesh_indices_to_bin(m, (/ i, j, k /)) + m % get_bin_from_indices([ i, j, k ]) score_index = 1 do l = 1, size(t % filter) diff --git a/src/cmfd_execute.F90 b/src/cmfd_execute.F90 index e87bff930c..b9aae0da5f 100644 --- a/src/cmfd_execute.F90 +++ b/src/cmfd_execute.F90 @@ -218,7 +218,7 @@ contains use error, only: warning, fatal_error use global, only: meshes, source_bank, work, n_user_meshes, cmfd use mesh_header, only: RegularMesh - use mesh, only: count_bank_sites, get_mesh_indices + use mesh, only: count_bank_sites use message_passing use string, only: to_str @@ -295,7 +295,7 @@ contains do i = 1, int(work,4) ! Determine spatial bin - call get_mesh_indices(m, source_bank(i) % xyz, ijk, in_mesh) + call m % get_indices(source_bank(i) % xyz, ijk, in_mesh) ! Determine energy bin n_groups = size(cmfd % egrid) - 1 diff --git a/src/mesh.F90 b/src/mesh.F90 index c25e3f7a6a..6d94ffba48 100644 --- a/src/mesh.F90 +++ b/src/mesh.F90 @@ -10,116 +10,6 @@ module mesh contains -!=============================================================================== -! GET_MESH_BIN determines the tally bin for a particle in a structured mesh -!=============================================================================== - - pure subroutine get_mesh_bin(m, xyz, bin) - type(RegularMesh), intent(in) :: m ! mesh pointer - real(8), intent(in) :: xyz(:) ! coordinates - integer, intent(out) :: bin ! tally bin - - integer :: n ! size of mesh - integer :: d ! mesh dimension index - integer :: ijk(3) ! indices in mesh - logical :: in_mesh ! was given coordinate in mesh at all? - - ! Get number of dimensions - n = m % n_dimension - - ! Loop over the dimensions of the mesh - do d = 1, n - - ! Check for cases where particle is outside of mesh - if (xyz(d) < m % lower_left(d)) then - bin = NO_BIN_FOUND - return - elseif (xyz(d) > m % upper_right(d)) then - bin = NO_BIN_FOUND - return - end if - end do - - ! Determine indices - call get_mesh_indices(m, xyz, ijk, in_mesh) - - ! Convert indices to bin - if (in_mesh) then - bin = mesh_indices_to_bin(m, ijk) - else - bin = NO_BIN_FOUND - end if - - end subroutine get_mesh_bin - -!=============================================================================== -! GET_MESH_INDICES determines the indices of a particle in a structured mesh -!=============================================================================== - - pure subroutine get_mesh_indices(m, xyz, ijk, in_mesh) - type(RegularMesh), intent(in) :: m - real(8), intent(in) :: xyz(:) ! coordinates to check - integer, intent(out) :: ijk(:) ! indices in mesh - logical, intent(out) :: in_mesh ! were given coords in mesh? - - ! Find particle in mesh - ijk(:m % n_dimension) = ceiling((xyz(:m % n_dimension) - m % lower_left)/m % width) - - ! Determine if particle is in mesh - if (any(ijk(:m % n_dimension) < 1) .or. & - any(ijk(:m % n_dimension) > m % dimension)) then - in_mesh = .false. - else - in_mesh = .true. - end if - - end subroutine get_mesh_indices - -!=============================================================================== -! MESH_INDICES_TO_BIN maps (i), (i,j), or (i,j,k) indices to a single bin number -! for use in a TallyObject results array -!=============================================================================== - - pure function mesh_indices_to_bin(m, ijk) result(bin) - type(RegularMesh), intent(in) :: m - integer, intent(in) :: ijk(:) - integer :: bin - - if (m % n_dimension == 1) then - bin = ijk(1) - elseif (m % n_dimension == 2) then - bin = (ijk(2) - 1) * m % dimension(1) + ijk(1) - elseif (m % n_dimension == 3) then - bin = ((ijk(3) - 1) * m % dimension(2) + (ijk(2) - 1)) & - * m % dimension(1) + ijk(1) - end if - - end function mesh_indices_to_bin - -!=============================================================================== -! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject results array to -! (i), (i,j), or (i,j,k) indices -!=============================================================================== - - pure subroutine bin_to_mesh_indices(m, bin, ijk) - type(RegularMesh), intent(in) :: m - integer, intent(in) :: bin - integer, intent(out) :: ijk(:) - - if (m % n_dimension == 1) then - ijk(1) = bin - else if (m % n_dimension == 2) then - ijk(1) = mod(bin - 1, m % dimension(1)) + 1 - ijk(2) = (bin - 1)/m % dimension(1) + 1 - else if (m % n_dimension == 3) then - ijk(1) = mod(bin - 1, m % dimension(1)) + 1 - ijk(2) = mod(bin - 1, m % dimension(1) * m % dimension(2)) & - / m % dimension(1) + 1 - ijk(3) = (bin - 1)/(m % dimension(1) * m % dimension(2)) + 1 - end if - - end subroutine bin_to_mesh_indices - !=============================================================================== ! COUNT_BANK_SITES determines the number of fission bank sites in each cell of a ! given mesh as well as an optional energy group structure. This can be used for @@ -172,7 +62,7 @@ contains ! loop over fission sites and count how many are in each mesh box FISSION_SITES: do i = 1, n_sites ! determine scoring bin for entropy mesh - call get_mesh_indices(m, bank_array(i) % xyz, ijk, in_mesh) + call m % get_indices(bank_array(i) % xyz, ijk, in_mesh) ! if outside mesh, skip particle if (.not. in_mesh) then @@ -215,226 +105,4 @@ contains end subroutine count_bank_sites -!=============================================================================== -! MESH_INTERSECTS determines if a line between xyz0 and xyz1 intersects the -! outer boundary of the given mesh. This is important for determining whether a -! track will score to a mesh tally. -!=============================================================================== - - pure function mesh_intersects_1d(m, xyz0, xyz1) result(intersects) - type(RegularMesh), intent(in) :: m - real(8), intent(in) :: xyz0(1) - real(8), intent(in) :: xyz1(1) - logical :: intersects - - real(8) :: x0 ! track start point - real(8) :: x1 ! track end point - real(8) :: xm0 ! lower-left coordinates of mesh - real(8) :: xm1 ! upper-right coordinates of mesh - - ! Copy coordinates of starting point - x0 = xyz0(1) - - ! Copy coordinates of ending point - x1 = xyz1(1) - - ! Copy coordinates of mesh lower_left - xm0 = m % lower_left(1) - - ! Copy coordinates of mesh upper_right - xm1 = m % upper_right(1) - - ! Set default value for intersects - intersects = .false. - - ! Check if line intersects left surface - if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then - intersects = .true. - return - end if - - ! Check if line intersects right surface - if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then - intersects = .true. - return - end if - - end function mesh_intersects_1d - - pure function mesh_intersects_2d(m, xyz0, xyz1) result(intersects) - type(RegularMesh), intent(in) :: m - real(8), intent(in) :: xyz0(2) - real(8), intent(in) :: xyz1(2) - logical :: intersects - - real(8) :: x0, y0 ! track start point - real(8) :: x1, y1 ! track end point - real(8) :: xi, yi ! track intersection point with mesh - real(8) :: xm0, ym0 ! lower-left coordinates of mesh - real(8) :: xm1, ym1 ! upper-right coordinates of mesh - - ! Copy coordinates of starting point - x0 = xyz0(1) - y0 = xyz0(2) - - ! Copy coordinates of ending point - x1 = xyz1(1) - y1 = xyz1(2) - - ! Copy coordinates of mesh lower_left - xm0 = m % lower_left(1) - ym0 = m % lower_left(2) - - ! Copy coordinates of mesh upper_right - xm1 = m % upper_right(1) - ym1 = m % upper_right(2) - - ! Set default value for intersects - intersects = .false. - - ! Check if line intersects left surface -- calculate the intersection point - ! y - if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then - yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0) - if (yi >= ym0 .and. yi < ym1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects back surface -- calculate the intersection point - ! x - if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then - xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0) - if (xi >= xm0 .and. xi < xm1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects right surface -- calculate the intersection - ! point y - if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then - yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0) - if (yi >= ym0 .and. yi < ym1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects front surface -- calculate the intersection point - ! x - if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then - xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0) - if (xi >= xm0 .and. xi < xm1) then - intersects = .true. - return - end if - end if - - end function mesh_intersects_2d - - pure function mesh_intersects_3d(m, xyz0, xyz1) result(intersects) - type(RegularMesh), intent(in) :: m - real(8), intent(in) :: xyz0(3) - real(8), intent(in) :: xyz1(3) - logical :: intersects - - real(8) :: x0, y0, z0 ! track start point - real(8) :: x1, y1, z1 ! track end point - real(8) :: xi, yi, zi ! track intersection point with mesh - real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh - real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh - - ! Copy coordinates of starting point - x0 = xyz0(1) - y0 = xyz0(2) - z0 = xyz0(3) - - ! Copy coordinates of ending point - x1 = xyz1(1) - y1 = xyz1(2) - z1 = xyz1(3) - - ! Copy coordinates of mesh lower_left - xm0 = m % lower_left(1) - ym0 = m % lower_left(2) - zm0 = m % lower_left(3) - - ! Copy coordinates of mesh upper_right - xm1 = m % upper_right(1) - ym1 = m % upper_right(2) - zm1 = m % upper_right(3) - - ! Set default value for intersects - intersects = .false. - - ! Check if line intersects left surface -- calculate the intersection point - ! (y,z) - if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then - yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0) - zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0) - if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects back surface -- calculate the intersection point - ! (x,z) - if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then - xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0) - zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0) - if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects bottom surface -- calculate the intersection - ! point (x,y) - if ((z0 < zm0 .and. z1 > zm0) .or. (z0 > zm0 .and. z1 < zm0)) then - xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0) - yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0) - if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects right surface -- calculate the intersection point - ! (y,z) - if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then - yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0) - zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0) - if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects front surface -- calculate the intersection point - ! (x,z) - if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then - xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0) - zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0) - if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then - intersects = .true. - return - end if - end if - - ! Check if line intersects top surface -- calculate the intersection point - ! (x,y) - if ((z0 < zm1 .and. z1 > zm1) .or. (z0 > zm1 .and. z1 < zm1)) then - xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0) - yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0) - if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then - intersects = .true. - return - end if - end if - - end function mesh_intersects_3d - end module mesh diff --git a/src/mesh_header.F90 b/src/mesh_header.F90 index 9da06f813e..9c40ccbef8 100644 --- a/src/mesh_header.F90 +++ b/src/mesh_header.F90 @@ -1,5 +1,7 @@ module mesh_header + use constants, only: NO_BIN_FOUND + implicit none !=============================================================================== @@ -16,6 +18,363 @@ module mesh_header real(8), allocatable :: lower_left(:) ! lower-left corner of mesh real(8), allocatable :: upper_right(:) ! upper-right corner of mesh real(8), allocatable :: width(:) ! width of each mesh cell + contains + procedure :: get_bin => regular_get_bin + procedure :: get_indices => regular_get_indices + procedure :: get_bin_from_indices => regular_get_bin_from_indices + procedure :: get_indices_from_bin => regular_get_indices_from_bin + procedure :: intersects => regular_intersects end type RegularMesh +contains + +!=============================================================================== +! GET_MESH_BIN determines the tally bin for a particle in a structured mesh +!=============================================================================== + + pure subroutine regular_get_bin(this, xyz, bin) + class(RegularMesh), intent(in) :: this + real(8), intent(in) :: xyz(:) ! coordinates + integer, intent(out) :: bin ! tally bin + + integer :: n ! size of mesh + integer :: d ! mesh dimension index + integer :: ijk(3) ! indices in mesh + logical :: in_mesh ! was given coordinate in mesh at all? + + ! Get number of dimensions + n = this % n_dimension + + ! Loop over the dimensions of the mesh + do d = 1, n + + ! Check for cases where particle is outside of mesh + if (xyz(d) < this % lower_left(d)) then + bin = NO_BIN_FOUND + return + elseif (xyz(d) > this % upper_right(d)) then + bin = NO_BIN_FOUND + return + end if + end do + + ! Determine indices + call this % get_indices(xyz, ijk, in_mesh) + + ! Convert indices to bin + if (in_mesh) then + bin = this % get_bin_from_indices(ijk) + else + bin = NO_BIN_FOUND + end if + + end subroutine regular_get_bin + +!=============================================================================== +! GET_MESH_INDICES determines the indices of a particle in a structured mesh +!=============================================================================== + + pure subroutine regular_get_indices(this, xyz, ijk, in_mesh) + class(RegularMesh), intent(in) :: this + real(8), intent(in) :: xyz(:) ! coordinates to check + integer, intent(out) :: ijk(:) ! indices in mesh + logical, intent(out) :: in_mesh ! were given coords in mesh? + + ! Find particle in mesh + ijk(:this % n_dimension) = ceiling((xyz(:this % n_dimension) - & + this % lower_left)/this % width) + + ! Determine if particle is in mesh + if (any(ijk(:this % n_dimension) < 1) .or. & + any(ijk(:this % n_dimension) > this % dimension)) then + in_mesh = .false. + else + in_mesh = .true. + end if + + end subroutine regular_get_indices + +!=============================================================================== +! MESH_INDICES_TO_BIN maps (i), (i,j), or (i,j,k) indices to a single bin number +! for use in a TallyObject results array +!=============================================================================== + + pure function regular_get_bin_from_indices(this, ijk) result(bin) + class(RegularMesh), intent(in) :: this + integer, intent(in) :: ijk(:) + integer :: bin + + if (this % n_dimension == 1) then + bin = ijk(1) + elseif (this % n_dimension == 2) then + bin = (ijk(2) - 1) * this % dimension(1) + ijk(1) + elseif (this % n_dimension == 3) then + bin = ((ijk(3) - 1) * this % dimension(2) + (ijk(2) - 1)) & + * this % dimension(1) + ijk(1) + end if + + end function regular_get_bin_from_indices + +!=============================================================================== +! BIN_TO_MESH_INDICES maps a single mesh bin from a TallyObject results array to +! (i), (i,j), or (i,j,k) indices +!=============================================================================== + + pure subroutine regular_get_indices_from_bin(this, bin, ijk) + class(RegularMesh), intent(in) :: this + integer, intent(in) :: bin + integer, intent(out) :: ijk(:) + + if (this % n_dimension == 1) then + ijk(1) = bin + else if (this % n_dimension == 2) then + ijk(1) = mod(bin - 1, this % dimension(1)) + 1 + ijk(2) = (bin - 1)/this % dimension(1) + 1 + else if (this % n_dimension == 3) then + ijk(1) = mod(bin - 1, this % dimension(1)) + 1 + ijk(2) = mod(bin - 1, this % dimension(1) * this % dimension(2)) & + / this % dimension(1) + 1 + ijk(3) = (bin - 1)/(this % dimension(1) * this % dimension(2)) + 1 + end if + + end subroutine regular_get_indices_from_bin + +!=============================================================================== +! MESH_INTERSECTS determines if a line between xyz0 and xyz1 intersects the +! outer boundary of the given mesh. This is important for determining whether a +! track will score to a mesh tally. +!=============================================================================== + + pure function regular_intersects(this, xyz0, xyz1) result(intersects) + class(RegularMesh), intent(in) :: this + real(8), intent(in) :: xyz0(1) + real(8), intent(in) :: xyz1(1) + logical :: intersects + + select case(this % n_dimension) + case (1) + intersects = mesh_intersects_1d(this, xyz0, xyz1) + case (2) + intersects = mesh_intersects_2d(this, xyz0, xyz1) + case (3) + intersects = mesh_intersects_3d(this, xyz0, xyz1) + end select + end function regular_intersects + + pure function mesh_intersects_1d(m, xyz0, xyz1) result(intersects) + type(RegularMesh), intent(in) :: m + real(8), intent(in) :: xyz0(:) + real(8), intent(in) :: xyz1(:) + logical :: intersects + + real(8) :: x0 ! track start point + real(8) :: x1 ! track end point + real(8) :: xm0 ! lower-left coordinates of mesh + real(8) :: xm1 ! upper-right coordinates of mesh + + ! Copy coordinates of starting point + x0 = xyz0(1) + + ! Copy coordinates of ending point + x1 = xyz1(1) + + ! Copy coordinates of mesh lower_left + xm0 = m % lower_left(1) + + ! Copy coordinates of mesh upper_right + xm1 = m % upper_right(1) + + ! Set default value for intersects + intersects = .false. + + ! Check if line intersects left surface + if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then + intersects = .true. + return + end if + + ! Check if line intersects right surface + if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then + intersects = .true. + return + end if + + end function mesh_intersects_1d + + pure function mesh_intersects_2d(m, xyz0, xyz1) result(intersects) + type(RegularMesh), intent(in) :: m + real(8), intent(in) :: xyz0(:) + real(8), intent(in) :: xyz1(:) + logical :: intersects + + real(8) :: x0, y0 ! track start point + real(8) :: x1, y1 ! track end point + real(8) :: xi, yi ! track intersection point with mesh + real(8) :: xm0, ym0 ! lower-left coordinates of mesh + real(8) :: xm1, ym1 ! upper-right coordinates of mesh + + ! Copy coordinates of starting point + x0 = xyz0(1) + y0 = xyz0(2) + + ! Copy coordinates of ending point + x1 = xyz1(1) + y1 = xyz1(2) + + ! Copy coordinates of mesh lower_left + xm0 = m % lower_left(1) + ym0 = m % lower_left(2) + + ! Copy coordinates of mesh upper_right + xm1 = m % upper_right(1) + ym1 = m % upper_right(2) + + ! Set default value for intersects + intersects = .false. + + ! Check if line intersects left surface -- calculate the intersection point + ! y + if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then + yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0) + if (yi >= ym0 .and. yi < ym1) then + intersects = .true. + return + end if + end if + + ! Check if line intersects back surface -- calculate the intersection point + ! x + if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then + xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0) + if (xi >= xm0 .and. xi < xm1) then + intersects = .true. + return + end if + end if + + ! Check if line intersects right surface -- calculate the intersection + ! point y + if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then + yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0) + if (yi >= ym0 .and. yi < ym1) then + intersects = .true. + return + end if + end if + + ! Check if line intersects front surface -- calculate the intersection point + ! x + if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then + xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0) + if (xi >= xm0 .and. xi < xm1) then + intersects = .true. + return + end if + end if + + end function mesh_intersects_2d + + pure function mesh_intersects_3d(m, xyz0, xyz1) result(intersects) + type(RegularMesh), intent(in) :: m + real(8), intent(in) :: xyz0(:) + real(8), intent(in) :: xyz1(:) + logical :: intersects + + real(8) :: x0, y0, z0 ! track start point + real(8) :: x1, y1, z1 ! track end point + real(8) :: xi, yi, zi ! track intersection point with mesh + real(8) :: xm0, ym0, zm0 ! lower-left coordinates of mesh + real(8) :: xm1, ym1, zm1 ! upper-right coordinates of mesh + + ! Copy coordinates of starting point + x0 = xyz0(1) + y0 = xyz0(2) + z0 = xyz0(3) + + ! Copy coordinates of ending point + x1 = xyz1(1) + y1 = xyz1(2) + z1 = xyz1(3) + + ! Copy coordinates of mesh lower_left + xm0 = m % lower_left(1) + ym0 = m % lower_left(2) + zm0 = m % lower_left(3) + + ! Copy coordinates of mesh upper_right + xm1 = m % upper_right(1) + ym1 = m % upper_right(2) + zm1 = m % upper_right(3) + + ! Set default value for intersects + intersects = .false. + + ! Check if line intersects left surface -- calculate the intersection point + ! (y,z) + if ((x0 < xm0 .and. x1 > xm0) .or. (x0 > xm0 .and. x1 < xm0)) then + yi = y0 + (xm0 - x0) * (y1 - y0) / (x1 - x0) + zi = z0 + (xm0 - x0) * (z1 - z0) / (x1 - x0) + if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then + intersects = .true. + return + end if + end if + + ! Check if line intersects back surface -- calculate the intersection point + ! (x,z) + if ((y0 < ym0 .and. y1 > ym0) .or. (y0 > ym0 .and. y1 < ym0)) then + xi = x0 + (ym0 - y0) * (x1 - x0) / (y1 - y0) + zi = z0 + (ym0 - y0) * (z1 - z0) / (y1 - y0) + if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then + intersects = .true. + return + end if + end if + + ! Check if line intersects bottom surface -- calculate the intersection + ! point (x,y) + if ((z0 < zm0 .and. z1 > zm0) .or. (z0 > zm0 .and. z1 < zm0)) then + xi = x0 + (zm0 - z0) * (x1 - x0) / (z1 - z0) + yi = y0 + (zm0 - z0) * (y1 - y0) / (z1 - z0) + if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then + intersects = .true. + return + end if + end if + + ! Check if line intersects right surface -- calculate the intersection point + ! (y,z) + if ((x0 < xm1 .and. x1 > xm1) .or. (x0 > xm1 .and. x1 < xm1)) then + yi = y0 + (xm1 - x0) * (y1 - y0) / (x1 - x0) + zi = z0 + (xm1 - x0) * (z1 - z0) / (x1 - x0) + if (yi >= ym0 .and. yi < ym1 .and. zi >= zm0 .and. zi < zm1) then + intersects = .true. + return + end if + end if + + ! Check if line intersects front surface -- calculate the intersection point + ! (x,z) + if ((y0 < ym1 .and. y1 > ym1) .or. (y0 > ym1 .and. y1 < ym1)) then + xi = x0 + (ym1 - y0) * (x1 - x0) / (y1 - y0) + zi = z0 + (ym1 - y0) * (z1 - z0) / (y1 - y0) + if (xi >= xm0 .and. xi < xm1 .and. zi >= zm0 .and. zi < zm1) then + intersects = .true. + return + end if + end if + + ! Check if line intersects top surface -- calculate the intersection point + ! (x,y) + if ((z0 < zm1 .and. z1 > zm1) .or. (z0 > zm1 .and. z1 < zm1)) then + xi = x0 + (zm1 - z0) * (x1 - x0) / (z1 - z0) + yi = y0 + (zm1 - z0) * (y1 - y0) / (z1 - z0) + if (xi >= xm0 .and. xi < xm1 .and. yi >= ym0 .and. yi < ym1) then + intersects = .true. + return + end if + end if + + end function mesh_intersects_3d + end module mesh_header diff --git a/src/output.F90 b/src/output.F90 index b4e620e134..4ea5810415 100644 --- a/src/output.F90 +++ b/src/output.F90 @@ -12,7 +12,6 @@ module output use global use math, only: t_percentile use mesh_header, only: RegularMesh - use mesh, only: mesh_indices_to_bin, bin_to_mesh_indices use message_passing, only: master, n_procs use nuclide_header use particle_header, only: LocalCoord, Particle @@ -1055,7 +1054,7 @@ contains do i = 1, n_cells ! Get the indices for this cell - call bin_to_mesh_indices(m, i, ijk) + call m % get_indices_from_bin(i, ijk) matches(i_filter_mesh) % bins % data(1) = i ! Write the header for this cell diff --git a/src/physics.F90 b/src/physics.F90 index f40b7fd7b8..8081aebcbd 100644 --- a/src/physics.F90 +++ b/src/physics.F90 @@ -8,7 +8,6 @@ module physics use global use material_header, only: Material use math - use mesh, only: get_mesh_indices use message_passing use nuclide_header use output, only: write_message @@ -1090,7 +1089,7 @@ contains if (ufs) then ! Determine indices on ufs mesh for current location - call get_mesh_indices(ufs_mesh, p % coord(1) % xyz, ijk, in_mesh) + call ufs_mesh % get_indices(p % coord(1) % xyz, ijk, in_mesh) if (.not. in_mesh) then call write_particle_restart(p) call fatal_error("Source site outside UFS mesh!") diff --git a/src/physics_mg.F90 b/src/physics_mg.F90 index b5d34838c2..c412e9e32a 100644 --- a/src/physics_mg.F90 +++ b/src/physics_mg.F90 @@ -8,7 +8,6 @@ module physics_mg use material_header, only: Material use math, only: rotate_angle use mgxs_header, only: Mgxs, MgxsContainer - use mesh, only: get_mesh_indices use message_passing use output, only: write_message use particle_header, only: Particle @@ -196,7 +195,7 @@ contains if (ufs) then ! Determine indices on ufs mesh for current location - call get_mesh_indices(ufs_mesh, p % coord(1) % xyz, ijk, in_mesh) + call ufs_mesh % get_indices(p % coord(1) % xyz, ijk, in_mesh) if (.not. in_mesh) then call write_particle_restart(p) diff --git a/src/plot.F90 b/src/plot.F90 index 067d0c56fc..b8ac3cfe62 100644 --- a/src/plot.F90 +++ b/src/plot.F90 @@ -10,8 +10,6 @@ module plot use geometry_header, only: Cell, root_universe use global use hdf5_interface - use mesh, only: get_mesh_indices - use mesh_header, only: RegularMesh use output, only: write_message, time_stamp use particle_header, only: LocalCoord, Particle use plot_header @@ -241,8 +239,8 @@ contains width = xyz_ur_plot - xyz_ll_plot associate (m => pl % meshlines_mesh) - call get_mesh_indices(m, xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh) - call get_mesh_indices(m, xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh) + call m % get_indices(xyz_ll_plot, ijk_ll(:m % n_dimension), in_mesh) + call m % get_indices(xyz_ur_plot, ijk_ur(:m % n_dimension), in_mesh) ! sweep through all meshbins on this plane and draw borders do i = ijk_ll(outer), ijk_ur(outer) diff --git a/src/tally.F90 b/src/tally.F90 index 338362c69f..cccf79d04a 100644 --- a/src/tally.F90 +++ b/src/tally.F90 @@ -9,10 +9,6 @@ module tally use geometry_header use global use math, only: t_percentile, calc_pn, calc_rn - use mesh, only: get_mesh_bin, bin_to_mesh_indices, & - get_mesh_indices, mesh_indices_to_bin, & - mesh_intersects_1d, mesh_intersects_2d, & - mesh_intersects_3d use mesh_header, only: RegularMesh use message_passing use output, only: header @@ -3259,19 +3255,13 @@ contains n_dim = m % n_dimension ! Determine indices for starting and ending location - call get_mesh_indices(m, xyz0, ijk0, start_in_mesh) - call get_mesh_indices(m, xyz1, ijk1, end_in_mesh) + call m % get_indices(xyz0, ijk0, start_in_mesh) + call m % get_indices(xyz1, ijk1, end_in_mesh) ! Check to see if start or end is in mesh -- if not, check if track still ! intersects with mesh if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then - if (n_dim == 1) then - if (.not. mesh_intersects_1d(m, xyz0, xyz1)) cycle - else if (n_dim == 2) then - if (.not. mesh_intersects_2d(m, xyz0, xyz1)) cycle - else - if (.not. mesh_intersects_3d(m, xyz0, xyz1)) cycle - end if + if (.not. m % intersects(xyz0, xyz1)) cycle end if ! Calculate number of surface crossings @@ -3350,7 +3340,7 @@ contains all(ijk0(:n_dim) <= m % dimension)) then filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 1 filter_matches(i_filter_mesh) % bins % data(1) = & - mesh_indices_to_bin(m, ijk0) + m % get_bin_from_indices(ijk0) filter_index = 1 do k = 1, size(t % filter) filter_index = filter_index + (filter_matches(t % & @@ -3389,7 +3379,7 @@ contains ijk0(d1) = ijk0(d1) + 1 filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 2 filter_matches(i_filter_mesh) % bins % data(1) = & - mesh_indices_to_bin(m, ijk0) + m % get_bin_from_indices(ijk0) filter_index = 1 do k = 1, size(t % filter) filter_index = filter_index + (filter_matches(t % & @@ -3412,7 +3402,7 @@ contains all(ijk0(:n_dim) <= m % dimension)) then filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 - 3 filter_matches(i_filter_mesh) % bins % data(1) = & - mesh_indices_to_bin(m, ijk0) + m % get_bin_from_indices(ijk0) filter_index = 1 do k = 1, size(t % filter) filter_index = filter_index + (filter_matches(t % & @@ -3451,7 +3441,7 @@ contains ijk0(d1) = ijk0(d1) - 1 filter_matches(i_filter_surf) % bins % data(1) = d1 * 4 filter_matches(i_filter_mesh) % bins % data(1) = & - mesh_indices_to_bin(m, ijk0) + m % get_bin_from_indices(ijk0) filter_index = 1 do k = 1, size(t % filter) filter_index = filter_index + (filter_matches(t % & diff --git a/src/tally_filter.F90 b/src/tally_filter.F90 index 0a7ce23cd4..f8283c0716 100644 --- a/src/tally_filter.F90 +++ b/src/tally_filter.F90 @@ -7,10 +7,6 @@ module tally_filter use global use hdf5_interface use mesh_header, only: RegularMesh - use mesh, only: get_mesh_bin, bin_to_mesh_indices, & - get_mesh_indices, mesh_indices_to_bin, & - mesh_intersects_1d, mesh_intersects_2d, & - mesh_intersects_3d use particle_header, only: Particle use string, only: to_str use tally_filter_header, only: TallyFilter, TallyFilterContainer, & @@ -266,7 +262,7 @@ contains if (estimator /= ESTIMATOR_TRACKLENGTH) then ! If this is an analog or collision tally, then there can only be one ! valid mesh bin. - call get_mesh_bin(m, p % coord(1) % xyz, bin) + call m % get_bin(p % coord(1) % xyz, bin) if (bin /= NO_BIN_FOUND) then call match % bins % push_back(bin) call match % weights % push_back(ONE) @@ -288,25 +284,13 @@ contains xyz1 = p % coord(1) % xyz - TINY_BIT * p % coord(1) % uvw ! Determine indices for starting and ending location. - call get_mesh_indices(m, xyz0, ijk0(:n), start_in_mesh) - call get_mesh_indices(m, xyz1, ijk1(:n), end_in_mesh) + call m % get_indices(xyz0, ijk0(:n), start_in_mesh) + call m % get_indices(xyz1, ijk1(:n), end_in_mesh) ! If this is the first iteration of the filter loop, check if the track ! intersects any part of the mesh. if ((.not. start_in_mesh) .and. (.not. end_in_mesh)) then - if (n == 1) then - if (.not. mesh_intersects_1d(m, xyz0, xyz1)) then - return - end if - else if (n == 2) then - if (.not. mesh_intersects_2d(m, xyz0, xyz1)) then - return - end if - else - if (.not. mesh_intersects_3d(m, xyz0, xyz1)) then - return - end if - end if + if (.not. m % intersects(xyz0, xyz1)) return end if ! ======================================================================== @@ -393,7 +377,7 @@ contains end if ! Assign the next tally bin and the score. - bin = mesh_indices_to_bin(m, ijk0(:n)) + bin = m % get_bin_from_indices(ijk0(:n)) call match % bins % push_back(bin) call match % weights % push_back(distance / total_distance) @@ -440,7 +424,7 @@ contains m => meshes(this % mesh) allocate(ijk(m % n_dimension)) - call bin_to_mesh_indices(m, bin, ijk) + call m % get_indices_from_bin(bin, ijk) if (m % n_dimension == 1) then label = "Mesh Index (" // trim(to_str(ijk(1))) // ")" elseif (m % n_dimension == 2) then diff --git a/src/trigger.F90 b/src/trigger.F90 index cb5f2c4c85..2927aea6c2 100644 --- a/src/trigger.F90 +++ b/src/trigger.F90 @@ -11,7 +11,6 @@ module trigger use global use string, only: to_str use output, only: warning, write_message - use mesh, only: bin_to_mesh_indices use mesh_header, only: RegularMesh use message_passing, only: master use trigger_header, only: TriggerObject @@ -344,7 +343,7 @@ contains do i = 1, n_cells ! Get the indices for this cell - call bin_to_mesh_indices(m, i, ijk) + call m % get_indices_from_bin(i, ijk) filter_matches(i_filter_mesh) % bins % data(1) = i do l = 1, n