adding mesh tally amalgamation algorithm

include/openmc/mesh.h

@@ -991,6 +991,9 @@ class LibMesh : public UnstructuredMesh {
 
   libMesh::MeshBase* mesh_ptr() const { return m_; };
 
+  //! setter for mesh tally amalgamtion
+  void set_mesh_tally_amalgamation(std::string cluster_element_integer_name);
+
 private:
   void initialize() override;
   void set_mesh_pointer_from_filename(const std::string& filename);
@@ -1030,6 +1033,19 @@ class LibMesh : public UnstructuredMesh {
                       //!< elements
   std::vector<int> elem_to_bin_map_; //!< mapping dof indices to bin indices for
                                      //!< active elements
+
+  bool amalgamation_ = false; //!< whether we are doing mesh and tally
+                              //!< amalgamation by default it's turned off.
+
+  int cluster_element_integer_index_ =
+    -1; //!< extra element integer index for element clustering
+  /*create a hash map where every element in a cluster would map to the first
+   * element of in that cluster if the element isn't part of a cluster then it
+   * will point to it self <any_element_in_a_cluster, first element in that
+   * cluster >
+   */
+  std::unordered_map<const libMesh::Elem*, const libMesh::Elem*>
+    clustering_element_mapping_;
 };
 
 #endif

src/mesh.cpp

@@ -3238,6 +3238,43 @@ LibMesh::LibMesh(libMesh::MeshBase& input_mesh, double length_multiplier)
   initialize();
 }
 
+void LibMesh::set_mesh_tally_amalgamation(
+  std::string cluster_element_integer_name)
+{
+
+  cluster_element_integer_index_ =
+    m_->has_elem_integer(cluster_element_integer_name)
+      ? m_->get_elem_integer_index(cluster_element_integer_name)
+      : -1;
+  amalgamation_ = (cluster_element_integer_index_ != -1);
+
+  // should we add a warning if amalgamation is false?
+
+  if (amalgamation_) {
+
+    // reseve the hash map for cluster elements
+    clustering_element_mapping_.reserve(m_->n_active_elem());
+
+    // adding clustering map
+    for (auto it = m_->active_elements_begin(); it != m_->active_elements_end();
+         it++) {
+
+      auto elem = *it;
+      auto cluster_elem = elem;
+      unsigned int cluster_id =
+        elem->get_extra_integer(cluster_element_integer_index_);
+
+      if (cluster_id != -1) {
+        auto first_element_in_a_cluster = m_->elem_ptr(cluster_id);
+
+        if (first_element_in_a_cluster and first_element_in_a_cluster->active())
+          cluster_elem = first_element_in_a_cluster;
+      }
+      clustering_element_mapping_.insert(std::make_pair(elem, cluster_elem));
+    }
+  }
+}
+
 // create the mesh from an input file
 LibMesh::LibMesh(const std::string& filename, double length_multiplier)
   : adaptive_(false)
@@ -3532,8 +3569,10 @@ int LibMesh::get_bin(Position r) const
 
 int LibMesh::get_bin_from_element(const libMesh::Elem* elem) const
 {
-  int bin =
-    adaptive_ ? elem_to_bin_map_[elem->id()] : elem->id() - first_element_id_;
+  auto tally_elem = amalgamation_ ? clustering_element_mapping_.at(elem) : elem;
+  int bin = adaptive_ ? elem_to_bin_map_[tally_elem->id()]
+                      : tally_elem->id() - first_element_id_;
+
   if (bin >= n_bins() || bin < 0) {
     fatal_error(fmt::format("Invalid bin: {}", bin));
   }