From efc77b74fe1ed7d70179581bdf7f9d4c3cbd0409 Mon Sep 17 00:00:00 2001 From: Sterling Harper Date: Fri, 16 Feb 2018 17:49:42 -0500 Subject: [PATCH] Move lattice get_indices and are_valid to C++ --- src/geometry_header.F90 | 177 ++++++++-------------------------------- src/lattice.cpp | 121 ++++++++++++++++++++++++++- src/lattice.h | 27 +++++- 3 files changed, 178 insertions(+), 147 deletions(-) diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index 0cc08f5e6..c58ab1d6c 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -95,6 +95,15 @@ module geometry_header integer(C_INT32_T) :: id end function lattice_id_c + function lattice_are_valid_indices_c(lat_ptr, i_xyz) & + bind(C, name='lattice_are_valid_indices') result (is_valid) + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: lat_ptr + integer(C_INT), intent(in) :: i_xyz(3) + logical(C_BOOL) :: is_valid + end function lattice_are_valid_indices_c + subroutine lattice_distance_c(lat_ptr, xyz, uvw, i_xyz, d, lattice_trans) & bind(C, name='lattice_distance') use ISO_C_BINDING @@ -107,6 +116,15 @@ module geometry_header integer(C_INT), intent(out) :: lattice_trans(3) end subroutine lattice_distance_c + subroutine lattice_get_indices_c(lat_ptr, xyz, i_xyz) & + bind(C, name='lattice_get_indices') + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: lat_ptr + real(C_DOUBLE), intent(in) :: xyz(3) + integer(C_INT), intent(out) :: i_xyz(3) + end subroutine lattice_get_indices_c + subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5') use ISO_C_BINDING use hdf5 @@ -146,37 +164,16 @@ module geometry_header contains procedure :: id => lattice_id + procedure :: are_valid_indices => lattice_are_valid_indices procedure :: distance => lattice_distance + procedure :: get_indices => lattice_get_indices procedure :: to_hdf5 => lattice_to_hdf5 - procedure(lattice_are_valid_indices_), deferred :: are_valid_indices - procedure(lattice_get_indices_), deferred :: get_indices procedure(lattice_get_local_xyz_), deferred :: get_local_xyz end type Lattice abstract interface -!=============================================================================== -! ARE_VALID_INDICES returns .true. if the given lattice indices fit within the -! bounds of the lattice. Returns false otherwise. - - function lattice_are_valid_indices_(this, i_xyz) result(is_valid) - import Lattice - class(Lattice), intent(in) :: this - integer, intent(in) :: i_xyz(3) - logical :: is_valid - end function lattice_are_valid_indices_ - -!=============================================================================== -! GET_INDICES returns the indices in a lattice for the given global xyz. - - function lattice_get_indices_(this, global_xyz) result(i_xyz) - import Lattice - class(Lattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer :: i_xyz(3) - end function lattice_get_indices_ - !=============================================================================== ! GET_LOCAL_XYZ returns the translated local version of the given global xyz. @@ -199,8 +196,6 @@ module geometry_header contains - procedure :: are_valid_indices => valid_inds_rect - procedure :: get_indices => get_inds_rect procedure :: get_local_xyz => get_local_rect end type RectLattice @@ -215,8 +210,6 @@ module geometry_header contains - procedure :: are_valid_indices => valid_inds_hex - procedure :: get_indices => get_inds_hex procedure :: get_local_xyz => get_local_hex end type HexLattice @@ -294,6 +287,13 @@ contains id = lattice_id_c(this % ptr) end function lattice_id + function lattice_are_valid_indices(this, i_xyz) result (is_valid) + class(Lattice), intent(in) :: this + integer(C_INT), intent(in) :: i_xyz(3) + logical(C_BOOL) :: is_valid + is_valid = lattice_are_valid_indices_c(this % ptr, i_xyz) + end function lattice_are_valid_indices + subroutine lattice_distance(this, xyz, uvw, i_xyz, d, lattice_trans) class(Lattice), intent(in) :: this real(C_DOUBLE), intent(in) :: xyz(3) @@ -304,128 +304,19 @@ contains call lattice_distance_c(this % ptr, xyz, uvw, i_xyz, d, lattice_trans) end subroutine lattice_distance + function lattice_get_indices(this, xyz) result(i_xyz) + class(Lattice), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3) + integer(C_INT) :: i_xyz(3) + call lattice_get_indices_c(this % ptr, xyz, i_xyz) + end function lattice_get_indices + subroutine lattice_to_hdf5(this, group) class(Lattice), intent(in) :: this integer(HID_T), intent(in) :: group call lattice_to_hdf5_c(this % ptr, group) end subroutine lattice_to_hdf5 -!=============================================================================== - - function valid_inds_rect(this, i_xyz) result(is_valid) - class(RectLattice), intent(in) :: this - integer, intent(in) :: i_xyz(3) - logical :: is_valid - - is_valid = all(i_xyz > 0 .and. i_xyz <= this % n_cells) - end function valid_inds_rect - -!=============================================================================== - - function valid_inds_hex(this, i_xyz) result(is_valid) - class(HexLattice), intent(in) :: this - integer, intent(in) :: i_xyz(3) - logical :: is_valid - - is_valid = (all(i_xyz > 0) .and. & - &i_xyz(1) < 2*this % n_rings .and. & - &i_xyz(2) < 2*this % n_rings .and. & - &i_xyz(1) + i_xyz(2) > this % n_rings .and. & - &i_xyz(1) + i_xyz(2) < 3*this % n_rings .and. & - &i_xyz(3) <= this % n_axial) - end function valid_inds_hex - -!=============================================================================== - - function get_inds_rect(this, global_xyz) result(i_xyz) - class(RectLattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer :: i_xyz(3) - - real(8) :: xyz(3) ! global_xyz alias - - xyz = global_xyz - - i_xyz(1) = ceiling((xyz(1) - this % lower_left(1))/this % pitch(1)) - i_xyz(2) = ceiling((xyz(2) - this % lower_left(2))/this % pitch(2)) - if (this % is_3d) then - i_xyz(3) = ceiling((xyz(3) - this % lower_left(3))/this % pitch(3)) - else - i_xyz(3) = 1 - end if - end function get_inds_rect - -!=============================================================================== - - function get_inds_hex(this, global_xyz) result(i_xyz) - class(HexLattice), intent(in) :: this - real(8), intent(in) :: global_xyz(3) - integer :: i_xyz(3) - - real(8) :: xyz(3) ! global xyz relative to the center - real(8) :: alpha ! Skewed coord axis - real(8) :: xyz_t(3) ! Local xyz - real(8) :: d, d_min ! Squared distance from cell centers - integer :: i, j, k ! Iterators - integer :: k_min ! Minimum distance index - - xyz(1) = global_xyz(1) - this % center(1) - xyz(2) = global_xyz(2) - this % center(2) - - ! Index z direction. - if (this % is_3d) then - xyz(3) = global_xyz(3) - this % center(3) - i_xyz(3) = ceiling(xyz(3)/this % pitch(2) + HALF*this % n_axial) - else - xyz(3) = global_xyz(3) - i_xyz(3) = 1 - end if - - ! Convert coordinates into skewed bases. The (x, alpha) basis is used to - ! find the index of the global coordinates to within 4 cells. - alpha = xyz(2) - xyz(1) / sqrt(THREE) - i_xyz(1) = floor(xyz(1) / (sqrt(THREE) / TWO * this % pitch(1))) - i_xyz(2) = floor(alpha / this % pitch(1)) - - ! Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but - ! the array is offset so that the indices never go below 1). - i_xyz(1) = i_xyz(1) + this % n_rings - i_xyz(2) = i_xyz(2) + this % n_rings - - ! Calculate the (squared) distance between the particle and the centers of - ! the four possible cells. Regular hexagonal tiles form a centroidal - ! Voronoi tessellation so the global xyz should be in the hexagonal cell - ! that it is closest to the center of. This method is used over a - ! method that uses the remainders of the floor divisions above because it - ! provides better finite precision performance. Squared distances are - ! used becasue they are more computationally efficient than normal - ! distances. - k = 1 - d_min = INFINITY - do i = 0, 1 - do j = 0, 1 - xyz_t = this % get_local_xyz(global_xyz, i_xyz + [j, i, 0]) - d = xyz_t(1)**2 + xyz_t(2)**2 - if (d < d_min) then - d_min = d - k_min = k - end if - k = k + 1 - end do - end do - - ! Select the minimum squared distance which corresponds to the cell the - ! coordinates are in. - if (k_min == 2) then - i_xyz(1) = i_xyz(1) + 1 - else if (k_min == 3) then - i_xyz(2) = i_xyz(2) + 1 - else if (k_min == 4) then - i_xyz(1) = i_xyz(1) + 1 - i_xyz(2) = i_xyz(2) + 1 - end if - end function get_inds_hex - !=============================================================================== function get_local_rect(this, global_xyz, i_xyz) result(local_xyz) diff --git a/src/lattice.cpp b/src/lattice.cpp index 473099034..82f15401d 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -131,6 +131,18 @@ RectLattice::RectLattice(pugi::xml_node lat_node) } } +//============================================================================== + +bool +RectLattice::are_valid_indices(const int i_xyz[3]) const +{ + return ( (i_xyz[0] > 0) && (i_xyz[0] <= n_cells[0]) + && (i_xyz[1] > 0) && (i_xyz[1] <= n_cells[1]) + && (i_xyz[2] > 0) && (i_xyz[2] <= n_cells[2])); +} + +//============================================================================== + std::pair> RectLattice::distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const @@ -191,6 +203,22 @@ RectLattice::distance(const double xyz[3], const double uvw[3], return {d, lattice_trans}; } +//============================================================================== + +std::array +RectLattice::get_indices(const double xyz[3]) const +{ + int ix {static_cast(std::ceil((xyz[0] - lower_left[0]) / pitch[0]))}; + int iy {static_cast(std::ceil((xyz[1] - lower_left[1]) / pitch[1]))}; + int iz; + if (is_3d) { + iz = static_cast(std::ceil((xyz[2] - lower_left[2]) / pitch[2])); + } else { + iz = 1; + } + return {ix, iy, iz}; +} + //============================================================================== // HexLattice implementation //============================================================================== @@ -229,7 +257,7 @@ HexLattice::HexLattice(pugi::xml_node lat_node) fatal_error("A hexagonal lattice with must have " "specified by 2 numbers."); } else if (!is_3d && (pitch_words.size() != 1)) { - fatal_error("A hexagonal lattice without must have
" + fatal_error("A hexagonal lattice without must have " "specified by 1 number."); } pitch[0] = stod(pitch_words[0]); @@ -342,6 +370,18 @@ HexLattice::HexLattice(pugi::xml_node lat_node) } } +//============================================================================== + +bool +HexLattice::are_valid_indices(const int i_xyz[3]) const +{ + return ((i_xyz[0] > 0) && (i_xyz[1] > 0) && (i_xyz[2] > 0) + && (i_xyz[0] < 2*n_rings) && (i_xyz[1] < 2*n_rings) + && (i_xyz[0] + i_xyz[1] > n_rings) + && (i_xyz[0] + i_xyz[1] < 3*n_rings) + && (i_xyz[2] <= n_axial)); +} + std::pair> HexLattice::distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const @@ -443,6 +483,74 @@ HexLattice::distance(const double xyz[3], const double uvw[3], return {d, lattice_trans}; } +//============================================================================== + +std::array +HexLattice::get_indices(const double xyz[3]) const +{ + // Offset the xyz by the lattice center. + double xyz_o[3] {xyz[0] - center[0], xyz[1] - center[1], xyz[2]}; + if (is_3d) {xyz_o[2] -= center[2];} + + // Index the z direction. + std::array out; + if (is_3d) { + out[2] = static_cast(std::ceil(xyz_o[2] / pitch[1] + 0.5 * n_axial)); + } else { + out[2] = 1; + } + + // Convert coordinates into skewed bases. The (x, alpha) basis is used to + // find the index of the global coordinates to within 4 cells. + double alpha = xyz_o[1] - xyz_o[0] / std::sqrt(3.0); + out[0] = static_cast(std::floor(xyz_o[0] + / (0.5*std::sqrt(3.0) * pitch[0]))); + out[1] = static_cast(std::floor(alpha / pitch[0])); + + // Add offset to indices (the center cell is (i_x, i_alpha) = (0, 0) but + // the array is offset so that the indices never go below 1). + out[0] += n_rings; + out[1] += n_rings; + + // Calculate the (squared) distance between the particle and the centers of + // the four possible cells. Regular hexagonal tiles form a Voronoi + // tessellation so the xyz should be in the hexagonal cell that it is closest + // to the center of. This method is used over a method that uses the + // remainders of the floor divisions above because it provides better finite + // precision performance. Squared distances are used becasue they are more + // computationally efficient than normal distances. + int k {1}; + int k_min {1}; + double d_min {INFTY}; + for (int i = 0; i < 2; i++) { + for (int j = 0; j < 2; j++) { + int i_xyz[3] {out[0] + j, out[1] + i, 0}; + std::array xyz_t = get_local_xyz(xyz, i_xyz); + double d = xyz_t[0]*xyz_t[0] + xyz_t[1]*xyz_t[1]; + if (d < d_min) { + d_min = d; + k_min = k; + } + k++; + } + } + + // Select the minimum squared distance which corresponds to the cell the + // coordinates are in. + if (k_min == 2) { + out[0] += 1; + } else if (k_min == 3) { + out[1] += 1; + } else if (k_min == 4) { + out[0] += 1; + out[1] += 1; + } + + return out; +} + +//============================================================================== + std::array HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const { @@ -500,6 +608,9 @@ extern "C" { int32_t lattice_id(Lattice *lat) {return lat->id;} + bool lattice_are_valid_indices(Lattice *lat, const int i_xyz[3]) + {return lat->are_valid_indices(i_xyz);} + void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3], const int i_xyz[3], double *d, int lattice_trans[3]) { @@ -510,6 +621,14 @@ extern "C" { lattice_trans[2] = ld.second[2]; } + void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3]) + { + std::array inds {lat->get_indices(xyz)}; + i_xyz[0] = inds[0]; + i_xyz[1] = inds[1]; + i_xyz[2] = inds[2]; + } + void lattice_to_hdf5(Lattice *lat, hid_t group) {lat->to_hdf5(group);} } diff --git a/src/lattice.h b/src/lattice.h index 2a4aa8865..81b0c2bdf 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -58,13 +58,26 @@ public: virtual ~Lattice() {} - //virtual bool are_valid_indices(const int i_xyz[3]) const = 0; + //! Check lattice indices. + //! @param i_xyz[3] The indices for a lattice tile. + //! @return true if the given indices fit within the lattice bounds. False + //! otherwise. + virtual bool are_valid_indices(const int i_xyz[3]) const = 0; + //! Find the next lattice surface crossing + //! @param xyz[3] A 3D Cartesian coordinate. + //! @param uvw[3] A 3D Cartesian direction. + //! @param i_xyz[3] The indices for a lattice tile. + //! @return The distance to the next crossing and an array indicating how the + //! lattice indices would change after crossing that boundary. virtual std::pair> distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const - = 0; + = 0; - //virtual void get_indices(const double global_xyz[3], int i_xyz[3]) const = 0; + //! Find the lattice tile indices for a given point. + //! @param xyz[3] A 3D Cartesian coordinate. + //! @return An array containing the indices of a lattice tile. + virtual std::array get_indices(const double xyz[3]) const = 0; //virtual void get_local_xyz(const double global_xyz[3], const int i_xyz[3], // double local_xyz[3]) const = 0; @@ -87,9 +100,13 @@ public: virtual ~RectLattice() {} + bool are_valid_indices(const int i_xyz[3]) const; + std::pair> distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const; + std::array get_indices(const double xyz[3]) const; + protected: std::array n_cells; //! Number of cells along each axis std::array lower_left; //! Global lower-left corner of the lattice @@ -103,9 +120,13 @@ public: virtual ~HexLattice() {} + bool are_valid_indices(const int i_xyz[3]) const; + std::pair> distance(const double xyz[3], const double uvw[3], const int i_xyz[3]) const; + std::array get_indices(const double xyz[3]) const; + std::array get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const;