diff --git a/src/geometry.F90 b/src/geometry.F90 index 1f241f66d4..a6965cffb9 100644 --- a/src/geometry.F90 +++ b/src/geometry.F90 @@ -422,70 +422,8 @@ contains LAT_TYPE: select type(lat) type is (RectLattice) - ! copy local coordinates - x = p % coord(j) % xyz(1) - y = p % coord(j) % xyz(2) - z = p % coord(j) % xyz(3) - - ! determine oncoming edge - x0 = sign(lat % pitch(1) * HALF, u) - y0 = sign(lat % pitch(2) * HALF, v) - - ! left and right sides - if (abs(x - x0) < FP_PRECISION) then - d = INFINITY - elseif (u == ZERO) then - d = INFINITY - else - d = (x0 - x)/u - end if - - d_lat = d - if (u > 0) then - level_lat_trans(:) = [1, 0, 0] - else - level_lat_trans(:) = [-1, 0, 0] - end if - - ! front and back sides - if (abs(y - y0) < FP_PRECISION) then - d = INFINITY - elseif (v == ZERO) then - d = INFINITY - else - d = (y0 - y)/v - end if - - if (d < d_lat) then - d_lat = d - if (v > 0) then - level_lat_trans(:) = [0, 1, 0] - else - level_lat_trans(:) = [0, -1, 0] - end if - end if - - if (lat % is_3d) then - z0 = sign(lat % pitch(3) * HALF, w) - - ! top and bottom sides - if (abs(z - z0) < FP_PRECISION) then - d = INFINITY - elseif (w == ZERO) then - d = INFINITY - else - d = (z0 - z)/w - end if - - if (d < d_lat) then - d_lat = d - if (w > 0) then - level_lat_trans(:) = [0, 0, 1] - else - level_lat_trans(:) = [0, 0, -1] - end if - end if - end if + call lat % distance(p % coord(j) % xyz, p % coord(j) % uvw, & + d_lat, level_lat_trans) type is (HexLattice) LAT_TYPE ! Copy local coordinates. diff --git a/src/geometry_header.F90 b/src/geometry_header.F90 index ad34c735ca..101bf8f9a1 100644 --- a/src/geometry_header.F90 +++ b/src/geometry_header.F90 @@ -95,6 +95,17 @@ module geometry_header integer(C_INT32_T) :: id end function lattice_id_c + subroutine lattice_distance_c(lat_ptr, xyz, uvw, d, lattice_trans) & + bind(C, name='lattice_distance') + use ISO_C_BINDING + implicit none + type(C_PTR), intent(in), value :: lat_ptr + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + real(C_DOUBLE), intent(out) :: d + integer(C_INT), intent(out) :: lattice_trans(3) + end subroutine lattice_distance_c + subroutine lattice_to_hdf5_c(lat_ptr, group) bind(C, name='lattice_to_hdf5') use ISO_C_BINDING use hdf5 @@ -134,6 +145,7 @@ module geometry_header contains procedure :: id => lattice_id + procedure :: distance => lattice_distance procedure :: to_hdf5 => lattice_to_hdf5 procedure(lattice_are_valid_indices_), deferred :: are_valid_indices @@ -281,6 +293,15 @@ contains id = lattice_id_c(this % ptr) end function lattice_id + subroutine lattice_distance(this, xyz, uvw, d, lattice_trans) + class(Lattice), intent(in) :: this + real(C_DOUBLE), intent(in) :: xyz(3) + real(C_DOUBLE), intent(in) :: uvw(3) + real(C_DOUBLE), intent(out) :: d + integer(C_INT), intent(out) :: lattice_trans(3) + call lattice_distance_c(this % ptr, xyz, uvw, d, lattice_trans) + end subroutine lattice_distance + subroutine lattice_to_hdf5(this, group) class(Lattice), intent(in) :: this integer(HID_T), intent(in) :: group diff --git a/src/lattice.cpp b/src/lattice.cpp index 006b7bbdd7..236e6d0bae 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -7,6 +7,9 @@ #include "hdf5_interface.h" #include "xml_interface.h" +//TODO: this is only inlcuded for constants that should be moved elsewhere +//#include "surface.h" + //TODO: remove this include #include @@ -81,9 +84,9 @@ RectLattice::RectLattice(pugi::xml_node lat_node) fatal_error("Number of entries on must be the same as the " "number of entries on ."); } - lower_left[0] = stoi(ll_words[0]); - lower_left[1] = stoi(ll_words[1]); - if (is_3d) {lower_left[2] = stoi(ll_words[2]);} + lower_left[0] = stod(ll_words[0]); + lower_left[1] = stod(ll_words[1]); + if (is_3d) {lower_left[2] = stod(ll_words[2]);} // Read the lattice pitches. std::string pitch_str{get_node_value(lat_node, "pitch")}; @@ -92,9 +95,9 @@ RectLattice::RectLattice(pugi::xml_node lat_node) fatal_error("Number of entries on must be the same as the " "number of entries on ."); } - pitch[0] = stoi(pitch_words[0]); - pitch[1] = stoi(pitch_words[1]); - if (is_3d) {pitch[2] = stoi(pitch_words[2]);} + pitch[0] = stod(pitch_words[0]); + pitch[1] = stod(pitch_words[1]); + if (is_3d) {pitch[2] = stod(pitch_words[2]);} // Read the universes and make sure the correct number was specified. int nx = n_cells[0]; @@ -129,6 +132,65 @@ RectLattice::RectLattice(pugi::xml_node lat_node) } } +std::pair> +RectLattice::distance(const double xyz[3], const double uvw[3]) const +{ + // Get short aliases to the coordinates. + double x {xyz[0]}; + double y {xyz[1]}; + double z {xyz[2]}; + double u {uvw[0]}; + double v {uvw[1]}; + + // Determine the oncoming edge. + double x0 {copysign(0.5 * pitch[0], u)}; + double y0 {copysign(0.5 * pitch[1], v)}; + + // Left and right sides + double d {INFTY}; + std::array lattice_trans; + if ((std::abs(x - x0) > FP_PRECISION) && u != 0) { + d = (x0 - x) / u; + if (u > 0) { + lattice_trans = {1, 0, 0}; + } else { + lattice_trans = {-1, 0, 0}; + } + } + + // Front and back sides + if ((std::abs(y - y0) > FP_PRECISION) && v != 0) { + double this_d = (y0 - y) / v; + if (this_d < d) { + d = this_d; + if (v > 0) { + lattice_trans = {0, 1, 0}; + } else { + lattice_trans = {0, -1, 0}; + } + } + } + + // Top and bottom sides + if (is_3d) { + double w {uvw[2]}; + double z0 {copysign(0.5 * pitch[2], w)}; + if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { + double this_d = (z0 - z) / w; + if (this_d < d) { + d = this_d; + if (w > 0) { + lattice_trans = {0, 0, 1}; + } else { + lattice_trans = {0, 0, -1}; + } + } + } + } + + return {d, lattice_trans}; +} + //============================================================================== // HexLattice implementation //============================================================================== @@ -138,6 +200,11 @@ HexLattice::HexLattice(pugi::xml_node lat_node) { } +std::pair> +HexLattice::distance(const double xyz[3], const double uvw[3]) const +{ +} + //============================================================================== extern "C" void @@ -173,6 +240,16 @@ extern "C" { int32_t lattice_id(Lattice *lat) {return lat->id;} + void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3], + double *d, int lattice_trans[3]) + { + std::pair> ld {lat->distance(xyz, uvw)}; + *d = ld.first; + lattice_trans[0] = ld.second[0]; + lattice_trans[1] = ld.second[1]; + lattice_trans[2] = ld.second[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 05b642e202..1825f396b2 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -3,6 +3,7 @@ #include #include +#include // For numeric_limits #include #include #include @@ -13,6 +14,14 @@ namespace openmc { +//============================================================================== +// Constants that should eventually be moved out of this file +//============================================================================== + +extern "C" double FP_PRECISION; +constexpr double INFTY{std::numeric_limits::max()}; +constexpr int C_NONE {-1}; + //============================================================================== // Constants //============================================================================== @@ -51,6 +60,9 @@ public: //virtual bool are_valid_indices(const int i_xyz[3]) const = 0; + virtual std::pair> + distance(const double xyz[3], const double uvw[3]) const = 0; + //virtual void get_indices(const double global_xyz[3], int i_xyz[3]) const = 0; //virtual void get_local_xyz(const double global_xyz[3], const int i_xyz[3], @@ -74,10 +86,13 @@ public: virtual ~RectLattice() {} + std::pair> + distance(const double xyz[3], const double uvw[3]) const; + protected: - std::array n_cells; //! Number of cells along each axis - std::array lower_left; //! Global lower-left corner of the lattice - std::array pitch; //! Lattice tile width along each axis + std::array n_cells; //! Number of cells along each axis + std::array lower_left; //! Global lower-left corner of the lattice + std::array pitch; //! Lattice tile width along each axis }; class HexLattice : public Lattice @@ -86,6 +101,9 @@ public: explicit HexLattice(pugi::xml_node lat_node); virtual ~HexLattice() {} + + std::pair> + distance(const double xyz[3], const double uvw[3]) const; }; } // namespace openmc