diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index e9519ca0a1..d2e3724972 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -71,18 +71,28 @@ public: //! Determine which bins were crossed by a particle // - //! \param[in] p Particle to check + //! \param[in] last_r Previous position of the particle + //! \param[in] r Current position of the particle + //! \param[in] u Particle direction //! \param[out] bins Bins that were crossed //! \param[out] lengths Fraction of tracklength in each bin - virtual void bins_crossed(const Particle& p, std::vector& bins, + virtual void bins_crossed(Position last_r, + Position r, + const Direction& u, + std::vector& bins, std::vector& lengths) const = 0; //! Determine which surface bins were crossed by a particle // - //! \param[in] p Particle to check + //! \param[in] r0 Previous position of the particle + //! \param[in] r1 Current position of the particle + //! \param[in] u Particle direction //! \param[out] bins Surface bins that were crossed virtual void - surface_bins_crossed(const Particle& p, std::vector& bins) const = 0; + surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const = 0; //! Get bin at a given position in space // @@ -137,7 +147,10 @@ public: int n_surface_bins() const override; - void bins_crossed(const Particle& p, std::vector& bins, + void bins_crossed(Position last_r, + Position r, + const Direction& u, + std::vector& bins, std::vector& lengths) const override; //! Count number of bank sites in each mesh bin / energy bin @@ -219,9 +232,11 @@ public: RegularMesh(pugi::xml_node node); // Overriden methods - - void surface_bins_crossed(const Particle& p, std::vector& bins) - const override; + void + surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const override; int get_index_in_direction(double r, int i) const override; @@ -259,9 +274,11 @@ public: RectilinearMesh(pugi::xml_node node); // Overriden methods - - void surface_bins_crossed(const Particle& p, std::vector& bins) - const override; + void + surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const override; int get_index_in_direction(double r, int i) const override; @@ -325,7 +342,7 @@ public: std::string bin_label(int bin) const override; void surface_bins_crossed(const Particle& p, - std::vector& bins) const override; + std::vector& bins) const; void to_hdf5(hid_t group) const override; @@ -348,10 +365,18 @@ public: MOABMesh(pugi::xml_node); MOABMesh(const std::string& filename); - void bins_crossed(const Particle& p, + void bins_crossed(Position last_r, + Position r, + const Direction& u, std::vector& bins, std::vector& lengths) const override; + void + surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const override; + int get_bin(Position r) const; int n_bins() const override; @@ -496,10 +521,18 @@ public: LibMesh(const std::string& filename); // Methods - void bins_crossed(const Particle& p, + void bins_crossed(Position last_r, + Position r, + const Direction& u, std::vector& bins, std::vector& lengths) const override; + void + surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const override; + int get_bin(Position r) const override; int n_bins() const override; diff --git a/include/openmc/tallies/filter_mesh.h b/include/openmc/tallies/filter_mesh.h index fa64877e36..8a712f4cb0 100644 --- a/include/openmc/tallies/filter_mesh.h +++ b/include/openmc/tallies/filter_mesh.h @@ -4,6 +4,7 @@ #include #include "openmc/tallies/filter.h" +#include "openmc/position.h" namespace openmc { @@ -42,11 +43,22 @@ public: virtual void set_mesh(int32_t mesh); + virtual void set_translation(const Position& translation); + + virtual void set_translation(const double translation[3]); + + virtual const Position& translation() const {return translation_;} + + virtual bool translated() const {return translated_;} + + protected: //---------------------------------------------------------------------------- // Data members int32_t mesh_; + bool translated_ {false}; + Position translation_ {0.0, 0.0, 0.0}; }; } // namespace openmc diff --git a/src/mesh.cpp b/src/mesh.cpp index 16ce630ae0..4274912f66 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -581,17 +581,15 @@ bool StructuredMesh::intersects_3d(Position& r0, Position r1, int* ijk) const return min_dist < INFTY; } -void StructuredMesh::bins_crossed(const Particle& p, std::vector& bins, - std::vector& lengths) const +void StructuredMesh::bins_crossed(Position last_r, + Position r, + const Direction& u, + std::vector& bins, + std::vector& lengths) const { // ======================================================================== // Determine where the track intersects the mesh and if it intersects at all. - // Copy the starting and ending coordinates of the particle. - Position last_r {p.r_last_}; - Position r {p.r()}; - Direction u {p.u()}; - // Compute the length of the entire track. double total_distance = (r - last_r).norm(); if (total_distance == 0.0) return; @@ -797,17 +795,12 @@ double RegularMesh::negative_grid_boundary(int* ijk, int i) const return lower_left_[i] + (ijk[i] - 1) * width_[i]; } -void RegularMesh::surface_bins_crossed(const Particle& p, - std::vector& bins) const +void +RegularMesh::surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const { - // ======================================================================== - // Determine if the track intersects the tally mesh. - - // Copy the starting and ending coordinates of the particle. - Position r0 {p.r_last_current_}; - Position r1 {p.r()}; - Direction u {p.u()}; - // Determine indices for starting and ending location. int n = n_dimension_; std::vector ijk0(n), ijk1(n); @@ -1107,17 +1100,11 @@ int RectilinearMesh::set_grid() return 0; } -void RectilinearMesh::surface_bins_crossed(const Particle& p, +void RectilinearMesh::surface_bins_crossed(Position r0, + Position r1, + const Direction& u, std::vector& bins) const { - // ======================================================================== - // Determine if the track intersects the tally mesh. - - // Copy the starting and ending coordinates of the particle. - Position r0 {p.r_last_current_}; - Position r1 {p.r()}; - Direction u {p.u()}; - // Determine indices for starting and ending location. int ijk0[3], ijk1[3]; bool start_in_mesh; @@ -1701,14 +1688,13 @@ MOABMesh::intersect_track(const moab::CartVect& start, } void -MOABMesh::bins_crossed(const Particle& p, - std::vector& bins, - std::vector& lengths) const +MOABUnstructuredMesh::bins_crossed(Position last_r, + Position r, + const Direction& u, + std::vector& bins, + std::vector& lengths) const { - Position last_r{p.r_last_}; - Position r{p.r()}; - Direction u{p.u()}; - u /= u.norm(); + // u /= u.norm(); moab::CartVect r0(last_r.x, last_r.y, last_r.z); moab::CartVect r1(r.x, r.y, r.z); moab::CartVect dir(u.x, u.y, u.z); @@ -1833,6 +1819,16 @@ double MOABMesh::tet_volume(moab::EntityHandle tet) const return 1.0 / 6.0 * (((p[1] - p[0]) * (p[2] - p[0])) % (p[3] - p[0])); } +void MOABUnstructuredMesh::surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const +{ + + // TODO: Implement triangle crossings here + throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."}; +} + int MOABMesh::get_bin(Position r) const { moab::EntityHandle tet = get_tet(r); @@ -2212,8 +2208,18 @@ int LibMesh::n_bins() const return m_->n_elem(); } -int LibMesh::n_surface_bins() const +void +LibMesh::surface_bins_crossed(Position r0, + Position r1, + const Direction& u, + std::vector& bins) const { + // TODO: Implement triangle crossings here + throw std::runtime_error{"Unstructured mesh surface tallies are not implemented."}; +} + +int +LibMesh::n_surface_bins() const { int n_bins = 0; for (int i = 0; i < this->n_bins(); i++) { const libMesh::Elem& e = get_element_from_bin(i); @@ -2300,7 +2306,9 @@ void LibMesh::write(const std::string& filename) const } void -LibMesh::bins_crossed(const Particle& p, +LibMesh::bins_crossed(Position last_r, + Position r, + const Direction& u, std::vector& bins, std::vector& lengths) const { diff --git a/src/tallies/filter_mesh.cpp b/src/tallies/filter_mesh.cpp index 7897173868..b85a08263e 100644 --- a/src/tallies/filter_mesh.cpp +++ b/src/tallies/filter_mesh.cpp @@ -34,14 +34,25 @@ void MeshFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { + + Position last_r = p.r_last_; + Position r = p.r(); + Position u = p.u(); + + // apply translation if present + if (translated_) { + last_r -= translation_; + r -= translation_; + } + if (estimator != TallyEstimator::TRACKLENGTH) { - auto bin = model::meshes[mesh_]->get_bin(p.r()); + auto bin = model::meshes[mesh_]->get_bin(r); if (bin >= 0) { match.bins_.push_back(bin); match.weights_.push_back(1.0); } } else { - model::meshes[mesh_]->bins_crossed(p, match.bins_, match.weights_); + model::meshes[mesh_]->bins_crossed(last_r, r, u, match.bins_, match.weights_); } } @@ -66,6 +77,18 @@ MeshFilter::set_mesh(int32_t mesh) n_bins_ = model::meshes[mesh_]->n_bins(); } +void MeshFilter::set_translation(const Position& translation) +{ + translated_ = true; + translation_ = translation; +} + +void MeshFilter::set_translation(const double translation[3]) +{ + translated_ = true; + translation_ = translation; +} + //============================================================================== // C-API functions //============================================================================== @@ -123,4 +146,27 @@ openmc_mesh_filter_set_mesh(int32_t index, int32_t index_mesh) return 0; } +extern "C" int +openmc_mesh_filter_set_translation(int32_t index, double translation[3]) +{ + // Make sure this is a valid index to an allocated filter + if (int err = verify_filter(index)) return err; + + // Get a pointer to the filter and downcast. + const auto& filt_base = model::tally_filters[index].get(); + auto* filt = dynamic_cast(filt_base); + + // Check the filter type + if (!filt) { + set_errmsg("Tried to set mesh on a non-mesh filter."); + return OPENMC_E_INVALID_TYPE; + } + + // Set the translation + filt->set_translation(translation); + + return 0; +} + + } // namespace openmc diff --git a/src/tallies/filter_meshsurface.cpp b/src/tallies/filter_meshsurface.cpp index 997fef7293..02fc815b75 100644 --- a/src/tallies/filter_meshsurface.cpp +++ b/src/tallies/filter_meshsurface.cpp @@ -11,7 +11,10 @@ void MeshSurfaceFilter::get_all_bins(const Particle& p, TallyEstimator estimator, FilterMatch& match) const { - model::meshes[mesh_]->surface_bins_crossed(p, match.bins_); + Position r0 = p.r_last_current_; + Position r1 = p.r(); + Direction u = p.u(); + model::meshes[mesh_]->surface_bins_crossed(r0, r1, u, match.bins_); for (auto b : match.bins_) match.weights_.push_back(1.0); }