diff --git a/include/openmc/mesh.h b/include/openmc/mesh.h
index 5c9272e93b9..f9060a277ad 100644
--- a/include/openmc/mesh.h
+++ b/include/openmc/mesh.h
@@ -1058,6 +1058,9 @@ class AdaptiveLibMesh : public LibMesh {
 
   void write(const std::string& filename) const override;
 
+  //! setter for mesh tally amalgamtion
+  void set_mesh_tally_amalgamation(std::string cluster_element_integer_name);
+
 protected:
   // Overridden methods
   int get_bin_from_element(const libMesh::Elem* elem) const override;
@@ -1073,6 +1076,21 @@ class AdaptiveLibMesh : public LibMesh {
                       //!< 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 clustering_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
diff --git a/src/mesh.cpp b/src/mesh.cpp
index 610d057cf08..d5a64ff9a9e 100644
--- a/src/mesh.cpp
+++ b/src/mesh.cpp
@@ -3741,6 +3741,44 @@ AdaptiveLibMesh::AdaptiveLibMesh(
   }
 }
 
+void AdaptiveLibMesh::set_mesh_tally_amalgamation(
+  std::string cluster_element_integer_name)
+{
+
+  clustering_element_integer_index_ =
+    m_->has_elem_integer(cluster_element_integer_name)
+      ? m_->get_elem_integer_index(cluster_element_integer_name)
+      : -1;
+  amalgamation_ = (clustering_element_integer_index_ != -1);
+
+  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(clustering_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));
+    }
+  } else {
+    fatal_error(fmt::format("No extra element integer named: {} found in the "
+                            "mesh",
+      cluster_element_integer_name));
+  }
+}
+
 int AdaptiveLibMesh::n_bins() const
 {
   return num_active_;
@@ -3770,7 +3808,9 @@ void AdaptiveLibMesh::write(const std::string& filename) const
 
 int AdaptiveLibMesh::get_bin_from_element(const libMesh::Elem* elem) const
 {
-  int bin = elem_to_bin_map_[elem->id()];
+  auto tally_elem = amalgamation_ ? clustering_element_mapping_.at(elem) : elem;
+  int bin = elem_to_bin_map_[tally_elem->id()];
+
   if (bin >= n_bins() || bin < 0) {
     fatal_error(fmt::format("Invalid bin: {}", bin));
   }