diff --git a/src/geometry.h b/src/geometry.h index 792d653ae..e069b9315 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -12,7 +12,7 @@ struct Position { Position() = default; Position(double x_, double y_, double z_) : x{x_}, y{y_}, z{z_} { }; - Position(double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; + Position(const double xyz[]) : x{xyz[0]}, y{xyz[1]}, z{xyz[2]} { }; Position& operator+=(Position); Position& operator+=(double); @@ -59,7 +59,7 @@ struct Angle : Position { double& w() { return z; } Angle() = default; Angle(double u, double v, double w) : Position{u, v, w} { }; - Angle(double uvw[]) : Position{uvw} { }; + Angle(const double uvw[]) : Position{uvw} { }; Angle(Position r) : Position{r} { }; }; diff --git a/src/hdf5_interface.h b/src/hdf5_interface.h index b0d2cc2c1..8b1f7a434 100644 --- a/src/hdf5_interface.h +++ b/src/hdf5_interface.h @@ -11,6 +11,8 @@ #include #include +#include "geometry.h" + namespace openmc { @@ -181,5 +183,12 @@ write_dataset(hid_t obj_id, const char* name, const std::array& buffer) write_dataset(obj_id, 1, dims, name, H5TypeMap::type_id, buffer.data(), false); } +template<> inline void +write_dataset(hid_t obj_id, const char* name, Position r) +{ + std::array buffer {r.x, r.y, r.z}; + write_dataset(obj_id, name, buffer); +} + } // namespace openmc #endif //HDF5_INTERFACE_H diff --git a/src/lattice.cpp b/src/lattice.cpp index c95762a63..c67fc419f 100644 --- a/src/lattice.cpp +++ b/src/lattice.cpp @@ -223,15 +223,14 @@ RectLattice::are_valid_indices(const int i_xyz[3]) const //============================================================================== std::pair> -RectLattice::distance(const double xyz[3], const double uvw[3], - const int i_xyz[3]) const +RectLattice::distance(Position r, Angle a, const int i_xyz[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]}; + double x = r.x; + double y = r.y; + double z = r.z; + double u = a.u(); + double v = a.v(); // Determine the oncoming edge. double x0 {copysign(0.5 * pitch[0], u)}; @@ -264,7 +263,7 @@ RectLattice::distance(const double xyz[3], const double uvw[3], // Top and bottom sides if (is_3d) { - double w {uvw[2]}; + double w {a.w()}; double z0 {copysign(0.5 * pitch[2], w)}; if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { double this_d = (z0 - z) / w; @@ -285,13 +284,13 @@ RectLattice::distance(const double xyz[3], const double uvw[3], //============================================================================== std::array -RectLattice::get_indices(const double xyz[3]) const +RectLattice::get_indices(Position r) const { - int ix {static_cast(std::ceil((xyz[0] - lower_left[0]) / pitch[0]))-1}; - int iy {static_cast(std::ceil((xyz[1] - lower_left[1]) / pitch[1]))-1}; + int ix {static_cast(std::ceil((r.x - lower_left.x) / pitch.x))-1}; + int iy {static_cast(std::ceil((r.y - lower_left.y) / pitch.y))-1}; int iz; if (is_3d) { - iz = static_cast(std::ceil((xyz[2] - lower_left[2]) / pitch[2]))-1; + iz = static_cast(std::ceil((r.z - lower_left.z) / pitch.z))-1; } else { iz = 0; } @@ -300,18 +299,15 @@ RectLattice::get_indices(const double xyz[3]) const //============================================================================== -std::array -RectLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const +Position +RectLattice::get_local_position(Position r, const int i_xyz[3]) const { - std::array local_xyz; - local_xyz[0] = global_xyz[0] - (lower_left[0] + (i_xyz[0] + 0.5)*pitch[0]); - local_xyz[1] = global_xyz[1] - (lower_left[1] + (i_xyz[1] + 0.5)*pitch[1]); + r.x -= (lower_left.x + (i_xyz[0] + 0.5)*pitch.x); + r.y -= (lower_left.y + (i_xyz[1] + 0.5)*pitch.y); if (is_3d) { - local_xyz[2] = global_xyz[2] - (lower_left[2] + (i_xyz[2] + 0.5)*pitch[2]); - } else { - local_xyz[2] = global_xyz[2]; + r.z -= (lower_left.z + (i_xyz[2] + 0.5)*pitch.z); } - return local_xyz; + return r; } //============================================================================== @@ -583,12 +579,11 @@ HexLattice::are_valid_indices(const int i_xyz[3]) const //============================================================================== std::pair> -HexLattice::distance(const double xyz[3], const double uvw[3], - const int i_xyz[3]) const +HexLattice::distance(Position r, Angle a, const int i_xyz[3]) const { // Compute the direction on the hexagonal basis. - double beta_dir = uvw[0] * std::sqrt(3.0) / 2.0 + uvw[1] / 2.0; - double gamma_dir = uvw[0] * std::sqrt(3.0) / 2.0 - uvw[1] / 2.0; + double beta_dir = a.u() * std::sqrt(3.0) / 2.0 + a.v() / 2.0; + double gamma_dir = a.u() * std::sqrt(3.0) / 2.0 - a.v() / 2.0; // Note that hexagonal lattice distance calculations are performed // using the particle's coordinates relative to the neighbor lattice @@ -600,15 +595,15 @@ HexLattice::distance(const double xyz[3], const double uvw[3], double d {INFTY}; std::array lattice_trans; double edge = -copysign(0.5*pitch[0], beta_dir); // Oncoming edge - std::array xyz_t; + Position r_t; if (beta_dir > 0) { const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1], i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } else { const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1], i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } - double beta = xyz_t[0] * std::sqrt(3.0) / 2.0 + xyz_t[1] / 2.0; + double beta = r_t.x * std::sqrt(3.0) / 2.0 + r_t.y / 2.0; if ((std::abs(beta - edge) > FP_PRECISION) && beta_dir != 0) { d = (edge - beta) / beta_dir; if (beta_dir > 0) { @@ -622,12 +617,12 @@ HexLattice::distance(const double xyz[3], const double uvw[3], edge = -copysign(0.5*pitch[0], gamma_dir); if (gamma_dir > 0) { const int i_xyz_t[3] {i_xyz[0]+1, i_xyz[1]-1, i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } else { const int i_xyz_t[3] {i_xyz[0]-1, i_xyz[1]+1, i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } - double gamma = xyz_t[0] * std::sqrt(3.0) / 2.0 - xyz_t[1] / 2.0; + double gamma = r_t.x * std::sqrt(3.0) / 2.0 - r_t.y / 2.0; if ((std::abs(gamma - edge) > FP_PRECISION) && gamma_dir != 0) { double this_d = (edge - gamma) / gamma_dir; if (this_d < d) { @@ -641,18 +636,18 @@ HexLattice::distance(const double xyz[3], const double uvw[3], } // Upper and lower sides. - edge = -copysign(0.5*pitch[0], uvw[1]); - if (uvw[1] > 0) { + edge = -copysign(0.5*pitch[0], a.v()); + if (a.v() > 0) { const int i_xyz_t[3] {i_xyz[0], i_xyz[1]+1, i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } else { const int i_xyz_t[3] {i_xyz[0], i_xyz[1]-1, i_xyz[2]}; - xyz_t = get_local_xyz(xyz, i_xyz_t); + r_t = get_local_position(r, i_xyz_t); } - if ((std::abs(xyz_t[1] - edge) > FP_PRECISION) && uvw[1] != 0) { - double this_d = (edge - xyz_t[1]) / uvw[1]; + if ((std::abs(r_t.y - edge) > FP_PRECISION) && a.v() != 0) { + double this_d = (edge - r_t.y) / a.v(); if (this_d < d) { - if (uvw[1] > 0) { + if (a.v() > 0) { lattice_trans = {0, 1, 0}; } else { lattice_trans = {0, -1, 0}; @@ -663,8 +658,8 @@ HexLattice::distance(const double xyz[3], const double uvw[3], // Top and bottom sides if (is_3d) { - double z {xyz[2]}; - double w {uvw[2]}; + double z = r.z; + double w = a.w(); double z0 {copysign(0.5 * pitch[1], w)}; if ((std::abs(z - z0) > FP_PRECISION) && w != 0) { double this_d = (z0 - z) / w; @@ -686,24 +681,24 @@ HexLattice::distance(const double xyz[3], const double uvw[3], //============================================================================== std::array -HexLattice::get_indices(const double xyz[3]) const +HexLattice::get_indices(Position r) 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];} + Position r_o {r.x - center.x, r.y - center.y, r.z}; + if (is_3d) {r_o.z -= center.z;} // 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))-1; + out[2] = static_cast(std::ceil(r_o.z / pitch[1] + 0.5 * n_axial))-1; } else { out[2] = 0; } // 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] + double alpha = r_o.y - r_o.x / std::sqrt(3.0); + out[0] = static_cast(std::floor(r_o.x / (0.5*std::sqrt(3.0) * pitch[0]))); out[1] = static_cast(std::floor(alpha / pitch[0])); @@ -725,8 +720,8 @@ HexLattice::get_indices(const double xyz[3]) const 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]; + Position r_t = get_local_position(r, i_xyz); + double d = r_t.x*r_t.x + r_t.y*r_t.y; if (d < d_min) { d_min = d; k_min = k; @@ -751,26 +746,19 @@ HexLattice::get_indices(const double xyz[3]) const //============================================================================== -std::array -HexLattice::get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const +Position +HexLattice::get_local_position(Position r, const int i_xyz[3]) const { - std::array local_xyz; - // x_l = x_g - (center + pitch_x*cos(30)*index_x) - local_xyz[0] = global_xyz[0] - (center[0] - + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]); + r.x -= (center.x + std::sqrt(3.0)/2.0 * (i_xyz[0] - n_rings + 1) * pitch[0]); // y_l = y_g - (center + pitch_x*index_x + pitch_y*sin(30)*index_y) - local_xyz[1] = global_xyz[1] - (center[1] - + (i_xyz[1] - n_rings + 1) * pitch[0] - + (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0); + r.y -= (center.y + (i_xyz[1] - n_rings + 1) * pitch[0] + + (i_xyz[0] - n_rings + 1) * pitch[0] / 2.0); if (is_3d) { - local_xyz[2] = global_xyz[2] - center[2] - + (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1]; - } else { - local_xyz[2] = global_xyz[2]; + r.z -= center.z - (0.5 * n_axial - i_xyz[2] - 0.5) * pitch[1]; } - return local_xyz; + return r; } //============================================================================== @@ -908,7 +896,9 @@ extern "C" { void lattice_distance(Lattice *lat, const double xyz[3], const double uvw[3], const int i_xyz[3], double *d, int lattice_trans[3]) { - std::pair> ld {lat->distance(xyz, uvw, i_xyz)}; + Position r {xyz}; + Angle a {uvw}; + std::pair> ld {lat->distance(r, a, i_xyz)}; *d = ld.first; lattice_trans[0] = ld.second[0]; lattice_trans[1] = ld.second[1]; @@ -917,7 +907,8 @@ extern "C" { void lattice_get_indices(Lattice *lat, const double xyz[3], int i_xyz[3]) { - std::array inds = lat->get_indices(xyz); + Position r {xyz}; + std::array inds = lat->get_indices(r); i_xyz[0] = inds[0]; i_xyz[1] = inds[1]; i_xyz[2] = inds[2]; @@ -926,10 +917,11 @@ extern "C" { void lattice_get_local_xyz(Lattice *lat, const double global_xyz[3], const int i_xyz[3], double local_xyz[3]) { - std::array xyz = lat->get_local_xyz(global_xyz, i_xyz); - local_xyz[0] = xyz[0]; - local_xyz[1] = xyz[1]; - local_xyz[2] = xyz[2]; + Position global {global_xyz}; + Position local = lat->get_local_position(global, i_xyz); + local_xyz[0] = local.x; + local_xyz[1] = local.y; + local_xyz[2] = local.z; } int32_t lattice_offset(Lattice *lat, int map, const int i_xyz[3]) diff --git a/src/lattice.h b/src/lattice.h index bc706a5a3..a6c449662 100644 --- a/src/lattice.h +++ b/src/lattice.h @@ -8,6 +8,7 @@ #include #include "constants.h" +#include "geometry.h" #include "hdf5.h" #include "pugixml.hpp" @@ -75,26 +76,26 @@ public: virtual bool are_valid_indices(const int i_xyz[3]) const = 0; //! \brief Find the next lattice surface crossing - //! @param xyz[3] A 3D Cartesian coordinate. - //! @param uvw[3] A 3D Cartesian direction. + //! @param r A 3D Cartesian coordinate. + //! @param a 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 + distance(Position r, Angle a, const int i_xyz[3]) const = 0; //! \brief Find the lattice tile indices for a given point. - //! @param xyz[3] A 3D Cartesian coordinate. + //! @param r 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 std::array get_indices(Position r) const = 0; //! \brief Get coordinates local to a lattice tile. - //! @param global_xyz[3] A 3D Cartesian coordinate. + //! @param r A 3D Cartesian coordinate. //! @param i_xyz[3] The indices for a lattice tile. //! @return Local 3D Cartesian coordinates. - virtual std::array - get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const = 0; + virtual Position + get_local_position(Position r, const int i_xyz[3]) const = 0; //! \brief Check flattened lattice index. //! @param indx The index for a lattice tile. @@ -193,12 +194,12 @@ public: 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; + distance(Position r, Angle a, const int i_xyz[3]) const; - std::array get_indices(const double xyz[3]) const; + std::array get_indices(Position r) const; - std::array - get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const; + Position + get_local_position(Position r, const int i_xyz[3]) const; int32_t& offset(int map, const int i_xyz[3]); @@ -207,9 +208,9 @@ public: void to_hdf5_inner(hid_t group_id) const; private: - 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 + Position lower_left; //!< Global lower-left corner of the lattice + Position pitch; //!< Lattice tile width along each axis // Convenience aliases int &nx {n_cells[0]}; @@ -233,12 +234,12 @@ public: 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; + distance(Position r, Angle a, const int i_xyz[3]) const; - std::array get_indices(const double xyz[3]) const; + std::array get_indices(Position r) const; - std::array - get_local_xyz(const double global_xyz[3], const int i_xyz[3]) const; + Position + get_local_position(Position r, const int i_xyz[3]) const; bool is_valid_index(int indx) const; @@ -251,7 +252,7 @@ public: private: int n_rings; //!< Number of radial tile positions int n_axial; //!< Number of axial tile positions - std::array center; //!< Global center of lattice + Position center; //!< Global center of lattice std::array pitch; //!< Lattice tile width and height };