diff --git a/include/openmc/geometry.h b/include/openmc/geometry.h index 9bb1a9ae7f..5feca2992a 100644 --- a/include/openmc/geometry.h +++ b/include/openmc/geometry.h @@ -1,6 +1,8 @@ #ifndef OPENMC_GEOMETRY_H #define OPENMC_GEOMETRY_H +#include +#include #include #include @@ -22,12 +24,22 @@ extern std::vector overlap_check_count; } // namespace model +//============================================================================== +// Information about nearest boundary crossing +//============================================================================== + +struct BoundaryInfo { + double distance {INFINITY}; //!< distance to nearest boundary + int surface_index {0}; //!< if boundary is surface, index in surfaces vector + int coord_level; //!< coordinate level after crossing boundary + std::array lattice_translation {}; //!< which way lattice indices will change +}; + //============================================================================== //! Check for overlapping cells at a particle's position. //============================================================================== -extern "C" bool -check_cell_overlap(Particle* p); +bool check_cell_overlap(Particle* p); //============================================================================== //! Locate a particle in the geometry tree and set its geometry data fields. @@ -41,23 +53,19 @@ check_cell_overlap(Particle* p); //! valid geometry coordinate stack. //============================================================================== -extern "C" bool -find_cell(Particle* p, bool use_neighbor_lists); +bool find_cell(Particle* p, bool use_neighbor_lists); //============================================================================== //! Move a particle into a new lattice tile. //============================================================================== -extern "C" void -cross_lattice(Particle* p, int lattice_translation[3]); +void cross_lattice(Particle* p, const BoundaryInfo& boundary); //============================================================================== //! Find the next boundary a particle will intersect. //============================================================================== -extern "C" void -distance_to_boundary(Particle* p, double* dist, int* surface_crossed, - int lattice_translation[3], int* next_level); +BoundaryInfo distance_to_boundary(Particle* p); } // namespace openmc diff --git a/src/geometry.cpp b/src/geometry.cpp index cde2cdc1b6..6eaa61fe24 100644 --- a/src/geometry.cpp +++ b/src/geometry.cpp @@ -32,8 +32,7 @@ std::vector overlap_check_count; // Non-member functions //============================================================================== -extern "C" bool -check_cell_overlap(Particle* p) +bool check_cell_overlap(Particle* p) { int n_coord = p->n_coord_; @@ -246,7 +245,7 @@ find_cell_inner(Particle* p, const NeighborList* neighbor_list) //============================================================================== -extern "C" bool +bool find_cell(Particle* p, bool use_neighbor_lists) { // Determine universe (if not yet set, use root universe). @@ -288,8 +287,8 @@ find_cell(Particle* p, bool use_neighbor_lists) //============================================================================== -extern "C" void -cross_lattice(Particle* p, int lattice_translation[3]) +void +cross_lattice(Particle* p, const BoundaryInfo& boundary) { auto& lat {*model::lattices[p->coord_[p->n_coord_-1].lattice]}; @@ -303,9 +302,9 @@ cross_lattice(Particle* p, int lattice_translation[3]) } // Set the lattice indices. - p->coord_[p->n_coord_-1].lattice_x += lattice_translation[0]; - p->coord_[p->n_coord_-1].lattice_y += lattice_translation[1]; - p->coord_[p->n_coord_-1].lattice_z += lattice_translation[2]; + p->coord_[p->n_coord_-1].lattice_x += boundary.lattice_translation[0]; + p->coord_[p->n_coord_-1].lattice_y += boundary.lattice_translation[1]; + p->coord_[p->n_coord_-1].lattice_z += boundary.lattice_translation[2]; std::array i_xyz {p->coord_[p->n_coord_-1].lattice_x, p->coord_[p->n_coord_-1].lattice_y, p->coord_[p->n_coord_-1].lattice_z}; @@ -346,16 +345,11 @@ cross_lattice(Particle* p, int lattice_translation[3]) //============================================================================== -extern "C" void -distance_to_boundary(Particle* p, double* dist, int* surface_crossed, - int lattice_translation[3], int* next_level) +BoundaryInfo distance_to_boundary(Particle* p) { - *dist = INFINITY; + BoundaryInfo info; double d_lat = INFINITY; double d_surf = INFINITY; - lattice_translation[0] = 0; - lattice_translation[1] = 0; - lattice_translation[2] = 0; int32_t level_surf_cross; std::array level_lat_trans; @@ -402,43 +396,43 @@ distance_to_boundary(Particle* p, double* dist, int* surface_crossed, // If the boundary on this coordinate level is coincident with a boundary on // a higher level then we need to make sure that the higher level boundary // is selected. This logic must consider floating point precision. + double& d = info.distance; if (d_surf < d_lat) { - if (*dist == INFINITY || ((*dist) - d_surf)/(*dist) >= FP_REL_PRECISION) { - *dist = d_surf; + if (d == INFINITY || (d - d_surf)/d >= FP_REL_PRECISION) { + d = d_surf; // If the cell is not simple, it is possible that both the negative and // positive half-space were given in the region specification. Thus, we // have to explicitly check which half-space the particle would be // traveling into if the surface is crossed if (c.simple_) { - *surface_crossed = level_surf_cross; + info.surface_index = level_surf_cross; } else { Position r_hit = r + d_surf * u; Surface& surf {*model::surfaces[std::abs(level_surf_cross)-1]}; Direction norm = surf.normal(r_hit); if (u.dot(norm) > 0) { - *surface_crossed = std::abs(level_surf_cross); + info.surface_index = std::abs(level_surf_cross); } else { - *surface_crossed = -std::abs(level_surf_cross); + info.surface_index = -std::abs(level_surf_cross); } } - lattice_translation[0] = 0; - lattice_translation[1] = 0; - lattice_translation[2] = 0; - *next_level = i + 1; + info.lattice_translation[0] = 0; + info.lattice_translation[1] = 0; + info.lattice_translation[2] = 0; + info.coord_level = i + 1; } } else { - if (*dist == INFINITY || ((*dist) - d_lat)/(*dist) >= FP_REL_PRECISION) { - *dist = d_lat; - *surface_crossed = F90_NONE; - lattice_translation[0] = level_lat_trans[0]; - lattice_translation[1] = level_lat_trans[1]; - lattice_translation[2] = level_lat_trans[2]; - *next_level = i + 1; + if (d == INFINITY || (d - d_lat)/d >= FP_REL_PRECISION) { + d = d_lat; + info.surface_index = 0; + info.lattice_translation = level_lat_trans; + info.coord_level = i + 1; } } } + return info; } //============================================================================== diff --git a/src/particle.cpp b/src/particle.cpp index b307aec8d1..4e0066b20e 100644 --- a/src/particle.cpp +++ b/src/particle.cpp @@ -202,12 +202,7 @@ Particle::transport() } // Find the distance to the nearest boundary - double d_boundary; - int surface_crossed; - int lattice_translation[3]; - int next_level; - distance_to_boundary(this, &d_boundary, &surface_crossed, - lattice_translation, &next_level); + auto boundary = distance_to_boundary(this); // Sample a distance to collision double d_collision; @@ -221,7 +216,7 @@ Particle::transport() } // Select smaller of the two distances - double distance = std::min(d_boundary, d_collision); + double distance = std::min(boundary.distance, d_collision); // Advance particle for (int j = 0; j < n_coord_; ++j) { @@ -244,11 +239,13 @@ Particle::transport() score_track_derivative(this, distance); } - if (d_collision > d_boundary) { + if (d_collision > boundary.distance) { // ==================================================================== // PARTICLE CROSSES SURFACE - if (next_level > 0) n_coord_ = next_level; + // Set surface that particle is on and adjust coordinate levels + surface_ = boundary.surface_index; + n_coord_ = boundary.coord_level; // Saving previous cell data for (int j = 0; j < n_coord_; ++j) { @@ -256,15 +253,14 @@ Particle::transport() } n_coord_last_ = n_coord_; - if (lattice_translation[0] != 0 || lattice_translation[1] != 0 || - lattice_translation[2] != 0) { + if (boundary.lattice_translation[0] != 0 || + boundary.lattice_translation[1] != 0 || + boundary.lattice_translation[2] != 0) { // Particle crosses lattice boundary - surface_ = 0; - cross_lattice(this, lattice_translation); + cross_lattice(this, boundary); event_ = EVENT_LATTICE; } else { // Particle crosses surface - surface_ = surface_crossed; this->cross_surface(); event_ = EVENT_SURFACE; }