Block restriction of libMesh unstructured mesh tallies (#3694)

Author: nuclearkevin

include/openmc/mesh.h

@@ -4,6 +4,7 @@
 #ifndef OPENMC_MESH_H
 #define OPENMC_MESH_H
 
+#include <set>
 #include <unordered_map>
 
 #include "hdf5.h"
@@ -969,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;
 
@@ -1007,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<unique_ptr<libMesh::PointLocatorBase>>
+    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();
@@ -1025,8 +1031,6 @@ class LibMesh : public UnstructuredMesh {
   unique_ptr<libMesh::MeshBase> unique_m_ =
     nullptr; //!< pointer to the libMesh MeshBase instance, only used if mesh is
              //!< created inside OpenMC
-  vector<unique_ptr<libMesh::PointLocatorBase>>
-    pl_; //!< per-thread point locators
   unique_ptr<libMesh::EquationSystems>
     equation_systems_; //!< pointer to the libMesh EquationSystems
                        //!< instance
@@ -1035,16 +1039,16 @@ class LibMesh : public UnstructuredMesh {
   std::unordered_map<std::string, unsigned int>
     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
 };
 
 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<libMesh::subdomain_id_type>& block_ids =
+      std::set<libMesh::subdomain_id_type>());
 
   // Overridden methods
   int n_bins() const override;
@@ -1056,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;
@@ -1064,6 +1070,9 @@ class AdaptiveLibMesh : public LibMesh {
 
 private:
   // Data members
+  const std::set<libMesh::subdomain_id_type>
+    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<libMesh::dof_id_type>

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<int> 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,22 +3732,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<libMesh::subdomain_id_type>& 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;
   }
@@ -3770,6 +3790,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()];