adding the base structure for on the fly amalgamation

include/openmc/mesh.h

@@ -991,6 +991,10 @@ 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 +1034,17 @@ 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)
@@ -3309,7 +3346,7 @@ void LibMesh::initialize()
     bin_to_elem_map_.reserve(m_->n_active_elem());
     elem_to_bin_map_.resize(m_->n_elem(), -1);
     for (auto it = m_->active_elements_begin(); it != m_->active_elements_end();
-         it++) {
+      it++) {
       auto elem = *it;
 
       bin_to_elem_map_.push_back(elem->id());
@@ -3325,6 +3362,7 @@ void LibMesh::initialize()
   upper_right_ = {ur(0), ur(1), ur(2)};
 }
 
+
 // Sample position within a tet for LibMesh type tets
 Position LibMesh::sample_element(int32_t bin, uint64_t* seed) const
 {
@@ -3468,7 +3506,7 @@ void LibMesh::set_score_data(const std::string& var_name,
   unsigned int std_dev_num = variable_map_.at(std_dev_name);
 
   for (auto it = m_->local_elements_begin(); it != m_->local_elements_end();
-       it++) {
+    it++) {
     if (!(*it)->active()) {
       continue;
     }
@@ -3532,8 +3570,9 @@ 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));
   }