diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h index 5c9272e93b9..e35a89b6897 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" @@ -971,7 +972,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; @@ -1009,16 +1010,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(); @@ -1027,8 +1033,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 @@ -1037,7 +1041,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 }; @@ -1045,8 +1048,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; @@ -1058,6 +1062,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; @@ -1066,6 +1072,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 c3a2d45818c..d1e15a2a553 100644 --- a/src/mesh.cpp +++ b/src/mesh.cpp @@ -3484,9 +3484,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_) { @@ -3536,7 +3533,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 @@ -3547,7 +3548,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 @@ -3688,6 +3693,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; @@ -3721,22 +3731,32 @@ 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, 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 elem = *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 (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; } @@ -3769,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()];