From ea8b727b8c3bfbe2db6936253bdb99cc34120319 Mon Sep 17 00:00:00 2001 From: GuySten Date: Mon, 31 Aug 2026 23:53:47 +0300 Subject: [PATCH] Let a mesh surface crossing change more than one index --- include/openmc/mesh.h | 47 +++++++++++++++++++--- src/mesh.cpp | 94 +++++++++++++++++++++++-------------------- 2 files changed, 92 insertions(+), 49 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 1c6044514bd..2c14bc74b2e 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -295,6 +295,13 @@ class Mesh { int n_dimension_ {-1}; //!< Number of dimensions }; +//! Maximum number of axes a structured mesh may have. +// +//! Mesh indices, shapes and the per-axis distance-to-surface arrays are all +//! sized by this, so it is the single place to change if a mesh type ever +//! needs more. +constexpr int MESH_MAX_AXES {3}; + class StructuredMesh : public Mesh { public: StructuredMesh() = default; @@ -302,14 +309,15 @@ class StructuredMesh : public Mesh { StructuredMesh(hid_t group) : Mesh {group} {}; virtual ~StructuredMesh() = default; - using MeshIndex = std::array; + using MeshIndex = std::array; struct MeshDistance { MeshDistance() = default; - MeshDistance(int _index, bool _max_surface, double _distance) - : next_index {_index}, max_surface {_max_surface}, distance {_distance} + MeshDistance(MeshIndex _offset, bool _max_surface, double _distance) + : offset {_offset}, max_surface {_max_surface}, distance {_distance} {} - int next_index {-1}; + //! Change in each index when this surface is crossed + MeshIndex offset {0, 0, 0}; bool max_surface {true}; double distance {INFTY}; bool operator<(const MeshDistance& o) const @@ -318,6 +326,12 @@ class StructuredMesh : public Mesh { } }; + //! Fold an index back into range on any periodic axis. + // + //! Called after an index is advanced by a MeshDistance offset. Meshes with a + //! periodic axis override this; for the rest it is a no-op. + virtual void sanitize_index(MeshIndex& idx) const {} + Position sample_element(int32_t bin, uint64_t* seed) const override { return sample_element(get_indices_from_bin(bin), seed); @@ -423,6 +437,17 @@ class StructuredMesh : public Mesh { virtual MeshDistance distance_to_grid_boundary(const MeshIndex& ijk, int i, const Position& r0, const Direction& u, double l) const = 0; + //! Distance to travel to re-enter the mesh from outside it + // + //! \param[in] ijk Current (out of range) indices + //! \param[in] distances Distance to the next surface on each axis + //! \param[in] traveled_distance Distance travelled so far + //! \param[out] k_max Axis whose surface is crossed last, or -1 + //! \return Updated travelled distance + double distance_to_mesh(const MeshIndex& ijk, + const std::array& distances, + double traveled_distance, int& k_max) const; + //! Get a label for the mesh bin std::string bin_label(int bin) const override; @@ -458,7 +483,8 @@ class StructuredMesh : public Mesh { virtual double volume(const MeshIndex& ijk) const = 0; // Data members - std::array shape_; //!< Number of mesh elements in each dimension + std::array + shape_; //!< Number of mesh elements in each dimension protected: }; @@ -587,6 +613,11 @@ class CylindricalMesh : public PeriodicStructuredMesh { CylindricalMesh(hid_t group); // Overridden methods + void sanitize_index(MeshIndex& idx) const override + { + idx[1] = sanitize_phi(idx[1]); + } + virtual MeshIndex get_indices(Position r, bool& in_mesh) const override; int get_index_in_direction(double r, int i) const override; @@ -654,6 +685,12 @@ class SphericalMesh : public PeriodicStructuredMesh { SphericalMesh(hid_t group); // Overridden methods + void sanitize_index(MeshIndex& idx) const override + { + idx[1] = sanitize_theta(idx[1]); + idx[2] = sanitize_phi(idx[2]); + } + virtual MeshIndex get_indices(Position r, bool& in_mesh) const override; int get_index_in_direction(double r, int i) const override; diff --git a/src/mesh.cpp b/src/mesh.cpp index 46698a6e390..f86649f8e1e 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -1208,7 +1208,7 @@ void StructuredMesh::raytrace_mesh( } // Calculate initial distances to next surfaces in all three dimensions - std::array distances; + std::array distances; for (int k = 0; k < n; ++k) { distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, 0.0); } @@ -1238,7 +1238,9 @@ void StructuredMesh::raytrace_mesh( // Update cell and calculate distance to next surface in k-direction. // The two other directions are still valid! - ijk[k] = distances[k].next_index; + for (int j = 0; j < MESH_MAX_AXES; ++j) + ijk[j] += distances[k].offset[j]; + sanitize_index(ijk); distances[k] = distance_to_grid_boundary(ijk, k, local_r, u, traveled_distance); @@ -1251,22 +1253,9 @@ void StructuredMesh::raytrace_mesh( tally.surface(ijk, k, !distances[k].max_surface, true); } else { // not inside mesh - - // For all directions outside the mesh, find the distance that we need - // to travel to reach the next surface. Use the largest distance, as - // only this will cross all outer surfaces. - int k_max {-1}; - for (int k = 0; k < n; ++k) { - if ((ijk[k] < 1 || ijk[k] > shape_[k]) && - (distances[k].distance > traveled_distance)) { - traveled_distance = distances[k].distance; - k_max = k; - } - } - // Assure some distance is traveled - if (k_max == -1) { - traveled_distance += TINY_BIT; - } + int k_max; + traveled_distance = + distance_to_mesh(ijk, distances, traveled_distance, k_max); // If r1 is not inside the mesh, exit here if (traveled_distance >= total_distance) @@ -1287,6 +1276,29 @@ void StructuredMesh::raytrace_mesh( } } +double StructuredMesh::distance_to_mesh(const MeshIndex& ijk, + const std::array& distances, + double traveled_distance, int& k_max) const +{ + // For all directions outside the mesh, find the distance that we need + // to travel to reach the next surface. Use the largest distance, as + // only this will cross all outer surfaces. + const int n = n_dimension_; + k_max = -1; + for (int k = 0; k < n; ++k) { + if ((ijk[k] < 1 || ijk[k] > shape_[k]) && + (distances[k].distance > traveled_distance)) { + traveled_distance = distances[k].distance; + k_max = k; + } + } + // Assure some distance is traveled + if (k_max == -1) { + traveled_distance += TINY_BIT; + } + return traveled_distance; +} + void StructuredMesh::bins_crossed(Position r0, Position r1, const Direction& u, vector& bins, vector& lengths) const { @@ -1527,16 +1539,15 @@ StructuredMesh::MeshDistance RegularMesh::distance_to_grid_boundary( double l) const { MeshDistance d; - d.next_index = ijk[i]; if (std::abs(u[i]) < FP_PRECISION) return d; d.max_surface = (u[i] > 0); if (d.max_surface && (ijk[i] <= shape_[i])) { - d.next_index++; + ++d.offset[i]; d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i]; } else if (!d.max_surface && (ijk[i] >= 1)) { - d.next_index--; + --d.offset[i]; d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i]; } @@ -1699,16 +1710,15 @@ StructuredMesh::MeshDistance RectilinearMesh::distance_to_grid_boundary( double l) const { MeshDistance d; - d.next_index = ijk[i]; if (std::abs(u[i]) < FP_PRECISION) return d; d.max_surface = (u[i] > 0); if (d.max_surface && (ijk[i] <= shape_[i])) { - d.next_index++; + ++d.offset[i]; d.distance = (positive_grid_boundary(ijk, i) - r0[i]) / u[i]; } else if (!d.max_surface && (ijk[i] > 0)) { - d.next_index--; + --d.offset[i]; d.distance = (negative_grid_boundary(ijk, i) - r0[i]) / u[i]; } return d; @@ -1966,18 +1976,16 @@ StructuredMesh::MeshDistance CylindricalMesh::find_z_crossing( const Position& r, const Direction& u, double l, int shell) const { MeshDistance d; - d.next_index = shell; - // Direction of flight is within xy-plane. Will never intersect z. if (std::abs(u.z) < FP_PRECISION) return d; d.max_surface = (u.z > 0.0); if (d.max_surface && (shell <= shape_[2])) { - d.next_index += 1; + ++d.offset[2]; d.distance = (grid_[2][shell] - r.z) / u.z; } else if (!d.max_surface && (shell > 0)) { - d.next_index -= 1; + --d.offset[2]; d.distance = (grid_[2][shell - 1] - r.z) / u.z; } return d; @@ -1990,15 +1998,14 @@ StructuredMesh::MeshDistance CylindricalMesh::distance_to_grid_boundary( if (i == 0) { return std::min( - MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])), - MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1))); + MeshDistance({1, 0, 0}, true, find_r_crossing(r0, u, l, ijk[i])), + MeshDistance({-1, 0, 0}, false, find_r_crossing(r0, u, l, ijk[i] - 1))); } else if (i == 1) { - return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true, - find_phi_crossing(r0, u, l, ijk[i])), - MeshDistance(sanitize_phi(ijk[i] - 1), false, - find_phi_crossing(r0, u, l, ijk[i] - 1))); + return std::min( + MeshDistance({0, 1, 0}, true, find_phi_crossing(r0, u, l, ijk[i])), + MeshDistance({0, -1, 0}, false, find_phi_crossing(r0, u, l, ijk[i] - 1))); } else { return find_z_crossing(r0, u, l, ijk[i]); @@ -2322,20 +2329,19 @@ StructuredMesh::MeshDistance SphericalMesh::distance_to_grid_boundary( if (i == 0) { return std::min( - MeshDistance(ijk[i] + 1, true, find_r_crossing(r0, u, l, ijk[i])), - MeshDistance(ijk[i] - 1, false, find_r_crossing(r0, u, l, ijk[i] - 1))); + MeshDistance({1, 0, 0}, true, find_r_crossing(r0, u, l, ijk[i])), + MeshDistance({-1, 0, 0}, false, find_r_crossing(r0, u, l, ijk[i] - 1))); } else if (i == 1) { - return std::min(MeshDistance(sanitize_theta(ijk[i] + 1), true, - find_theta_crossing(r0, u, l, ijk[i])), - MeshDistance(sanitize_theta(ijk[i] - 1), false, - find_theta_crossing(r0, u, l, ijk[i] - 1))); + return std::min( + MeshDistance({0, 1, 0}, true, find_theta_crossing(r0, u, l, ijk[i])), + MeshDistance( + {0, -1, 0}, false, find_theta_crossing(r0, u, l, ijk[i] - 1))); } else { - return std::min(MeshDistance(sanitize_phi(ijk[i] + 1), true, - find_phi_crossing(r0, u, l, ijk[i])), - MeshDistance(sanitize_phi(ijk[i] - 1), false, - find_phi_crossing(r0, u, l, ijk[i] - 1))); + return std::min( + MeshDistance({0, 0, 1}, true, find_phi_crossing(r0, u, l, ijk[i])), + MeshDistance({0, 0, -1}, false, find_phi_crossing(r0, u, l, ijk[i] - 1))); } }