From 2a555e5e9c417be4a7d88db4fecef2cee2b8fa65 Mon Sep 17 00:00:00 2001 From: nuclearkevin <66632997+nuclearkevin@users.noreply.github.com> Date: Fri, 19 Dec 2025 13:07:08 -0600 Subject: [PATCH 1/3] Block restriction. --- include/openmc/mesh.h | 9 +++++++-- src/mesh.cpp | 21 ++++++++++++++++----- 2 files changed, 23 insertions(+), 7 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 7ceb623dae1..3aa418a9238 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -4,6 +4,7 @@ #ifndef OPENMC_MESH_H #define OPENMC_MESH_H +#include #include #include "hdf5.h" @@ -1043,8 +1044,9 @@ class LibMesh : public UnstructuredMesh { class AdaptiveLibMesh : public LibMesh { public: // Constructor - AdaptiveLibMesh( - libMesh::MeshBase& input_mesh, double length_multiplier = 1.0); + AdaptiveLibMesh(libMesh::MeshBase& input_mesh, double length_multiplier = 1.0, + const std::set& block_ids = + std::set()); // Overridden methods int n_bins() const override; @@ -1064,6 +1066,9 @@ class AdaptiveLibMesh : public LibMesh { private: // Data members + const std::set + block_ids_; //!< subdomains of the mesh to tally on + const bool block_restrict_; //!< whether a subset of the mesh is being used const libMesh::dof_id_type num_active_; //!< cached number of active elements std::vector diff --git a/src/mesh.cpp b/src/mesh.cpp index 9a7a8e7517b..a434464fe21 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -3725,17 +3725,28 @@ double LibMesh::volume(int bin) const return this->get_element_from_bin(bin).volume(); } -AdaptiveLibMesh::AdaptiveLibMesh( - libMesh::MeshBase& input_mesh, double length_multiplier) - : LibMesh(input_mesh, length_multiplier), num_active_(m_->n_active_elem()) +AdaptiveLibMesh::AdaptiveLibMesh(libMesh::MeshBase& input_mesh, + double length_multiplier, + const std::set& block_ids) + : LibMesh(input_mesh, length_multiplier), block_ids_(block_ids), + block_restrict_(!block_ids_.empty()), + num_active_( + block_restrict_ + ? std::distance(m_->active_subdomain_set_elements_begin(block_ids_), + m_->active_subdomain_set_elements_end(block_ids_)) + : m_->n_active_elem()) { // if the mesh is adaptive elements aren't guaranteed by libMesh to be // contiguous in ID space, so we need to map from bin indices (defined over // active elements) to global dof ids bin_to_elem_map_.reserve(num_active_); elem_to_bin_map_.resize(m_->n_elem(), -1); - for (auto it = m_->active_elements_begin(); it != m_->active_elements_end(); - it++) { + auto begin = block_restrict_ + ? m_->active_subdomain_set_elements_begin(block_ids_) + : m_->active_elements_begin(); + auto end = block_restrict_ ? m_->active_subdomain_set_elements_end(block_ids_) + : m_->active_elements_end(); + for (auto it = begin; it != end; it++) { auto elem = *it; bin_to_elem_map_.push_back(elem->id()); From 6d1f305af5fce0fb63f9f248e8f9c28f2a58e826 Mon Sep 17 00:00:00 2001 From: nuclearkevin <66632997+nuclearkevin@users.noreply.github.com> Date: Fri, 19 Dec 2025 17:19:41 -0600 Subject: [PATCH 2/3] Add block restriction to point locators and improve scaling. --- include/openmc/mesh.h | 14 ++++++++----- src/mesh.cpp | 47 ++++++++++++++++++++++++++++++++++--------- 2 files changed, 47 insertions(+), 14 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 3aa418a9238..1681c06877c 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -970,7 +970,7 @@ class LibMesh : public UnstructuredMesh { Position sample_element(int32_t bin, uint64_t* seed) const override; - int get_bin(Position r) const override; + virtual int get_bin(Position r) const override; int n_bins() const override; @@ -1008,16 +1008,21 @@ class LibMesh : public UnstructuredMesh { protected: // Methods - //! Translate a bin value to an element reference virtual const libMesh::Elem& get_element_from_bin(int bin) const; //! Translate an element pointer to a bin index virtual int get_bin_from_element(const libMesh::Elem* elem) const; + // Data members libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set //!< during intialization + vector> + pl_; //!< per-thread point locators + libMesh::BoundingBox bbox_; //!< bounding box of the mesh + private: + // Methods void initialize() override; void set_mesh_pointer_from_filename(const std::string& filename); void build_eqn_sys(); @@ -1026,8 +1031,6 @@ class LibMesh : public UnstructuredMesh { unique_ptr unique_m_ = nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is //!< created inside OpenMC - vector> - pl_; //!< per-thread point locators unique_ptr equation_systems_; //!< pointer to the libMesh EquationSystems //!< instance @@ -1036,7 +1039,6 @@ class LibMesh : public UnstructuredMesh { std::unordered_map variable_map_; //!< mapping of variable names (tally scores) to libMesh //!< variable numbers - libMesh::BoundingBox bbox_; //!< bounding box of the mesh libMesh::dof_id_type first_element_id_; //!< id of the first element in the mesh }; @@ -1058,6 +1060,8 @@ class AdaptiveLibMesh : public LibMesh { void write(const std::string& filename) const override; + int get_bin(Position r) const override; + protected: // Overridden methods int get_bin_from_element(const libMesh::Elem* elem) const override; diff --git a/src/mesh.cpp b/src/mesh.cpp index a434464fe21..2d67217f755 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -3485,9 +3485,6 @@ void LibMesh::initialize() // assuming that unstructured meshes used in OpenMC are 3D n_dimension_ = 3; - if (length_multiplier_ > 0.0) { - libMesh::MeshTools::Modification::scale(*m_, length_multiplier_); - } // if OpenMC is managing the libMesh::MeshBase instance, prepare the mesh. // Otherwise assume that it is prepared by its owning application if (unique_m_) { @@ -3537,7 +3534,11 @@ Position LibMesh::centroid(int bin) const { const auto& elem = this->get_element_from_bin(bin); auto centroid = elem.vertex_average(); - return {centroid(0), centroid(1), centroid(2)}; + if (length_multiplier_ > 0.0) { + return length_multiplier_ * Position(centroid(0), centroid(1), centroid(2)); + } else { + return {centroid(0), centroid(1), centroid(2)}; + } } int LibMesh::n_vertices() const @@ -3548,7 +3549,11 @@ int LibMesh::n_vertices() const Position LibMesh::vertex(int vertex_id) const { const auto node_ref = m_->node_ref(vertex_id); - return {node_ref(0), node_ref(1), node_ref(2)}; + if (length_multiplier_ > 0.0) { + return length_multiplier_ * Position(node_ref(0), node_ref(1), node_ref(2)); + } else { + return {node_ref(0), node_ref(1), node_ref(2)}; + } } std::vector LibMesh::connectivity(int elem_id) const @@ -3689,6 +3694,11 @@ int LibMesh::get_bin(Position r) const // look-up a tet using the point locator libMesh::Point p(r.x, r.y, r.z); + if (length_multiplier_ > 0.0) { + // Scale the point down + p /= length_multiplier_; + } + // quick rejection check if (!bbox_.contains_point(p)) { return -1; @@ -3722,7 +3732,7 @@ const libMesh::Elem& LibMesh::get_element_from_bin(int bin) const double LibMesh::volume(int bin) const { - return this->get_element_from_bin(bin).volume(); + return this->get_element_from_bin(bin).volume() * length_multiplier_ * length_multiplier_ * length_multiplier_; } AdaptiveLibMesh::AdaptiveLibMesh(libMesh::MeshBase& input_mesh, @@ -3746,9 +3756,7 @@ AdaptiveLibMesh::AdaptiveLibMesh(libMesh::MeshBase& input_mesh, : m_->active_elements_begin(); auto end = block_restrict_ ? m_->active_subdomain_set_elements_end(block_ids_) : m_->active_elements_end(); - for (auto it = begin; it != end; it++) { - auto elem = *it; - + for (const auto & elem : libMesh::as_range(begin, end)) { bin_to_elem_map_.push_back(elem->id()); elem_to_bin_map_[elem->id()] = bin_to_elem_map_.size() - 1; } @@ -3781,6 +3789,27 @@ void AdaptiveLibMesh::write(const std::string& filename) const this->id_)); } +int AdaptiveLibMesh::get_bin(Position r) const +{ + // look-up a tet using the point locator + libMesh::Point p(r.x, r.y, r.z); + + if (length_multiplier_ > 0.0) { + // Scale the point down + p /= length_multiplier_; + } + + // quick rejection check + if (!bbox_.contains_point(p)) { + return -1; + } + + const auto& point_locator = pl_.at(thread_num()); + + const auto elem_ptr = (*point_locator)(p, &block_ids_); + return elem_ptr ? get_bin_from_element(elem_ptr) : -1; +} + int AdaptiveLibMesh::get_bin_from_element(const libMesh::Elem* elem) const { int bin = elem_to_bin_map_[elem->id()]; From e81ad7fce9d76e84cf4f81ecdf3080c01ad9e015 Mon Sep 17 00:00:00 2001 From: nuclearkevin <66632997+nuclearkevin@users.noreply.github.com> Date: Fri, 16 Jan 2026 13:36:00 -0600 Subject: [PATCH 3/3] Style changes. --- include/openmc/mesh.h | 2 +- src/mesh.cpp | 5 +++-- 2 files changed, 4 insertions(+), 3 deletions(-) diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 1681c06877c..2936d3ab1bd 100644 --- a/include/openmc/mesh.h +++ b/include/openmc/mesh.h @@ -1018,7 +1018,7 @@ class LibMesh : public UnstructuredMesh { libMesh::MeshBase* m_; //!< pointer to libMesh MeshBase instance, always set //!< during intialization vector> - pl_; //!< per-thread point locators + pl_; //!< per-thread point locators libMesh::BoundingBox bbox_; //!< bounding box of the mesh private: diff --git a/src/mesh.cpp b/src/mesh.cpp index 2d67217f755..6ee15b9c18b 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -3732,7 +3732,8 @@ const libMesh::Elem& LibMesh::get_element_from_bin(int bin) const double LibMesh::volume(int bin) const { - return this->get_element_from_bin(bin).volume() * length_multiplier_ * length_multiplier_ * length_multiplier_; + return this->get_element_from_bin(bin).volume() * length_multiplier_ * + length_multiplier_ * length_multiplier_; } AdaptiveLibMesh::AdaptiveLibMesh(libMesh::MeshBase& input_mesh, @@ -3756,7 +3757,7 @@ AdaptiveLibMesh::AdaptiveLibMesh(libMesh::MeshBase& input_mesh, : m_->active_elements_begin(); auto end = block_restrict_ ? m_->active_subdomain_set_elements_end(block_ids_) : m_->active_elements_end(); - for (const auto & elem : libMesh::as_range(begin, end)) { + for (const auto& elem : libMesh::as_range(begin, end)) { bin_to_elem_map_.push_back(elem->id()); elem_to_bin_map_[elem->id()] = bin_to_elem_map_.size() - 1; }