SF_mesh_utils_emi.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------
 
#if WITH_EMI_MODEL
 
#ifndef _SF_MESH_UTILS_EMI_H
#define _SF_MESH_UTILS_EMI_H
 
#include <algorithm>
#include <cstdio>
#include <cstdlib>
 
#include "SF_mesh_utils.h"
#include "SF_mesh_io_emi.h"
 
namespace SF {
 
 
const int MAX_INTERSECTIONS = 20;
// Adaptive merge-rank settings for EMI mesh merging.
// For communicator sizes up to the node-aware cap, we preserve the original
// all-ranks behavior. Once size > merge_rank_cap, we reduce the number of
// merge ranks based on the estimated global communication volume, bounded by
// that cap, and spread those merge ranks across the communicator. The lower
// bound avoids concentrating all merge work on too few ranks. These values are
// heuristics and may need retuning on different machines.
const int EMI_MIN_MERGE_RANKS = 8;
const int EMI_MAX_MERGE_RANKS_PER_NODE = 2;
const unsigned long long EMI_TARGET_BYTES_PER_MERGE_RANK = 8ull * 1024ull * 1024ull;
 
inline int emi_node_aware_merge_rank_cap(MPI_Comm comm)
{
  int size = 0;
  MPI_Comm_size(comm, &size);
 
  // Split the communicator into shared-memory subcommunicators, one per node.
  MPI_Comm node_comm = MPI_COMM_NULL;
  MPI_Comm_split_type(comm, MPI_COMM_TYPE_SHARED, 0, MPI_INFO_NULL, &node_comm);
 
  // Rank 0 in each node subcommunicator acts as the node representative.
  int node_rank = 0;
  MPI_Comm_rank(node_comm, &node_rank);
 
  int local_leader = (node_rank == 0) ? 1 : 0;
  int num_nodes = 0;
  // Count how many node representatives exist globally, which gives the number of nodes.
  MPI_Allreduce(&local_leader, &num_nodes, 1, MPI_INT, MPI_SUM, comm);
 
  MPI_Comm_free(&node_comm);
 
  // Allow only a small number of merge ranks per node, bounded by the communicator size.
  return std::min(size, std::max(EMI_MIN_MERGE_RANKS, num_nodes * EMI_MAX_MERGE_RANKS_PER_NODE));
}
 
template<class K, class TokenFn>
inline int emi_select_merge_destinations(const vector<K>& key_vec,
                                         size_t dsize,
                                         size_t bytes_per_entry,
                                         vector<int>& dest,
                                         MPI_Comm comm,
                                         const char* label,
                                         TokenFn token_fn)
{
  int size = 0, rank = 0;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  // Estimate how many bytes this helper exchanges locally and globally.
  const unsigned long long local_bytes =
      static_cast<unsigned long long>(dsize) * static_cast<unsigned long long>(bytes_per_entry);
  unsigned long long global_bytes = 0;
  MPI_Allreduce(&local_bytes, &global_bytes, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, comm);
 
  // Compute the node-aware upper bound on how many merge ranks we allow.
  const int merge_rank_cap = emi_node_aware_merge_rank_cap(comm);
  int nmerge = size;
  if (size > merge_rank_cap) {
    // For large runs, estimate how many merge ranks are needed from the total data volume.
    const int nmerge_floor = EMI_MIN_MERGE_RANKS;
    const int nmerge_cap = merge_rank_cap;
    const unsigned long long nmerge_est_ull =
        std::max(1ull, (global_bytes + EMI_TARGET_BYTES_PER_MERGE_RANK - 1ull) /
                           EMI_TARGET_BYTES_PER_MERGE_RANK);
    const int nmerge_est = static_cast<int>(std::min<unsigned long long>(nmerge_cap, nmerge_est_ull));
    // Clamp the estimate to the allowed range for this communicator.
    nmerge = std::max(nmerge_floor, nmerge_est);
  }
 
  dest.resize(dsize);
  for (size_t i = 0; i < dsize; ++i) {
    if (nmerge == size) {
      // Preserve the original all-ranks merge behavior for small and medium runs.
      dest[i] = static_cast<int>(token_fn(key_vec[i]) % static_cast<size_t>(size));
    } else {
      // First pick one logical merge bucket in [0, nmerge).
      const size_t bucket = token_fn(key_vec[i]) % static_cast<size_t>(nmerge);
      // Then spread those buckets over the full communicator instead of packing them
      // into ranks [0, nmerge), which reduces concentration on low ranks.
      dest[i] = static_cast<int>((bucket * static_cast<size_t>(size)) / static_cast<size_t>(nmerge));
    }
  }
 
  if (rank == 0 && std::getenv("OPENCARP_EMI_LOG_MERGE_RANKS") != NULL) {
    const char* mode = (nmerge == size) ? "all-ranks" : "limited";
    std::fprintf(stdout,
                 "[EMI merge ranks] %s: size=%d dsize=%zu global_bytes=%llu cap=%d nmerge=%d mode=%s\n",
                 label, size, dsize, global_bytes, merge_rank_cap, nmerge, mode);
    std::fflush(stdout);
  }
 
  return nmerge;
}
 
struct intersection_tags {                                                                                                       
  int tags[MAX_INTERSECTIONS];
  intersection_tags() {
    std::fill_n(tags,  MAX_INTERSECTIONS, -1);
  }
 
};
 
struct intersection_data {
  int tags[MAX_INTERSECTIONS];
  int data[MAX_INTERSECTIONS];
  int ranks[MAX_INTERSECTIONS];
 
  intersection_data() {
      std::fill_n(tags,  MAX_INTERSECTIONS, -1);
      std::fill_n(data,  MAX_INTERSECTIONS, -1);
      std::fill_n(ranks, MAX_INTERSECTIONS, -1);
  }
};
 
struct intersection_indices {
  mesh_int_t indices[MAX_INTERSECTIONS];
  intersection_indices() {
    std::fill_n(indices,  MAX_INTERSECTIONS, -1);
  }
};
template<class K, class V> inline
void assign_counter_face(hashmap::unordered_map<K, V>  & map, const MPI_Comm comm)
{
  size_t dsize = map.size();
  vector<K>    key_vec   (dsize);
  vector<V>    value_vec (dsize);
 
  // make key and value vector of elements without counterparts
  // which are from different ranks
  size_t idx=0;
  for(const auto & v : map) {
    if (v.second.second.eidx == -1){
      key_vec[idx]   = v.first;
      value_vec[idx] = v.second;
      idx++;
    }
  }
  dsize = idx;
  key_vec.resize(dsize);
  value_vec.resize(dsize);
 
  vector<int> perm, dest;
  emi_select_merge_destinations(key_vec, dsize, sizeof(K) + sizeof(V), dest, comm,
                                "assign_counter_face",
                                [](const K& key) { return hashmap::hash_ops<K>::hash(key); });
 
  commgraph<size_t> grph;
  grph.configure(dest, comm);
  size_t nrecv = sum(grph.rcnt);
 
  interval(perm, 0, dsize);
  binary_sort_copy(dest, perm);
 
  // fill send buffer and communicate
  vector<K> sbuff_key(dsize),   rbuff_key(nrecv);
  vector<V> sbuff_value(dsize), ibuff_value(dsize), rbuff_value(nrecv);
  for(size_t i=0; i<dsize; i++) {
    sbuff_key[i]   = key_vec[perm[i]];
    sbuff_value[i] = value_vec[perm[i]];
  }
  MPI_Exchange(grph, sbuff_key, rbuff_key, comm);
  MPI_Exchange(grph, sbuff_value, rbuff_value, comm);
 
  // Merge values for identical keys on the receiving side.
  // Each merged entry should contain both sides and preserve owner index_unique.
  hashmap::unordered_map<K,V>  rmap;
 
  // Merge all received entries for the same key.
  for (size_t i = 0; i < nrecv; i++) {
    auto it = rmap.find(rbuff_key[i]);
    if(it != rmap.end()) {
      // Add counter faces only if the tag numbers of the two faces are different.
      // That is, the counter face must have a different tag number from the current face.
      // However, in some special cases, both faces may belong to the extracellular domain,
      // which are removed later. 
      if (it->second.first.tag != rbuff_value[i].first.tag) {
        it->second.second    = rbuff_value[i].first;
      }
    } else {
      rmap.insert({rbuff_key[i], rbuff_value[i]});
    }
  }
 
  for (size_t i = 0; i < nrecv; i++ ) {
    auto it = rmap.find(rbuff_key[i]);
    if(it != rmap.end()) rbuff_value[i] = it->second;
  }
 
  // Send merged values back to the originating ranks.
  grph.transpose();
  MPI_Exchange(grph, rbuff_value, ibuff_value, comm);
 
  for (size_t i = 0; i < dsize; i++ ) {
    auto it = map.find(sbuff_key[i]);
    if (it != map.end()) it->second = ibuff_value[i];
  }
}
 
template<class K, class V> inline
void assign_unique_first_face(hashmap::unordered_map<K, V>  & map, const MPI_Comm comm)
{
  // Exchange and merge interface-face pairs across ranks.
  // Each key represents one physical interface; after this call, both sides
  // (first/second) carry rank/index data, and the owner side has index_unique set.
  //
  // Steps:
  // 1) Collect only partial entries (missing counter-face) that need MPI exchange.
  // 2) Hash-partition these keys to a merge owner rank.
  // 3) Exchange keys/values to owners and merge both sides for each key.
  // 4) Send merged entries back so every rank has complete data for its keys.
  size_t dsize = map.size();
  vector<K>    key_vec   (dsize);
  vector<V>    value_vec (dsize);
 
  // Collect only entries without counterparts (cross-rank interfaces).
  size_t idx=0;
  for(const auto & v : map) {
    if (v.second.second.index_unique == -1 && v.second.second.index_one == -1){
      key_vec[idx]   = v.first;
      value_vec[idx] = v.second;
      idx++;
    }
  }
  dsize = idx;
  key_vec.resize(dsize);
  value_vec.resize(dsize);
 
  vector<int> perm, dest;
  emi_select_merge_destinations(key_vec, dsize, sizeof(K) + sizeof(V), dest, comm,
                                "assign_unique_first_face",
                                [](const K& key) { return hashmap::hash_ops<K>::hash(key); });
 
  commgraph<size_t> grph;
  grph.configure(dest, comm);
  size_t nrecv = sum(grph.rcnt);
 
  interval(perm, 0, dsize);
  binary_sort_copy(dest, perm);
 
  // Send keys and partial values to the owner rank for each key.
  vector<K> sbuff_key(dsize),   rbuff_key(nrecv);
  vector<V> sbuff_value(dsize), ibuff_value(dsize), rbuff_value(nrecv);
  for(size_t i=0; i<dsize; i++) {
    sbuff_key[i]   = key_vec[perm[i]];
    sbuff_value[i] = value_vec[perm[i]];
  }
  MPI_Exchange(grph, sbuff_key, rbuff_key, comm);
  MPI_Exchange(grph, sbuff_value, rbuff_value, comm);
 
  // rmap accumulates merged entries per key on the receiving rank.
  // Multiple ranks may send the same key; we combine them into one pair.
  hashmap::unordered_map<K,V>  rmap;
 
  // Merge all received entries for the same key.
  for (size_t i = 0; i < nrecv; i++) {
    auto it = rmap.find(rbuff_key[i]);
    if(it != rmap.end()) {
      // Merge entries without losing existing owner info.
      auto& existing = it->second;
      using Face = std::decay_t<decltype(existing.first)>;
 
      auto update_slot = [&](Face& slot, const Face& src) {
        // Merge strategy for one side:
        // - index_unique is only set by the owner and must never be overwritten by a non-owner.
        // - rank/index_one are filled lazily when missing.
        if (src.index_unique >= 0) {
          slot.index_unique = src.index_unique;
          slot.rank = src.rank;
        } else {
          if (slot.rank < 0) slot.rank = src.rank;
        }
        if (slot.index_one < 0 && src.index_one >= 0) slot.index_one = src.index_one;
      };
 
      auto same_face = [&](const Face& a, const Face& b) {
        // Use (rank,index_one) to identify a side uniquely (index_one is local).
        return (a.rank >= 0 && b.rank >= 0 &&
                a.index_one >= 0 && b.index_one >= 0 &&
                a.rank == b.rank && a.index_one == b.index_one);
      };
 
      auto ingest = [&](const Face& incoming) {
        // Ignore empty placeholders (no useful info).
        if (incoming.index_one < 0 && incoming.index_unique < 0 && incoming.rank < 0) return;
        // Try to match incoming with an existing side (same rank or same index_one).
        if (same_face(existing.first, incoming) || existing.first.rank == incoming.rank) {
          update_slot(existing.first, incoming);
        } else if (same_face(existing.second, incoming) || existing.second.rank == incoming.rank) {
          update_slot(existing.second, incoming);
        } else if (existing.second.index_one == -1 && existing.second.index_unique == -1) {
          // If the second slot is empty, place incoming there (new side).
          existing.second = incoming;
        } else if (existing.first.index_one == -1 && existing.first.index_unique == -1) {
          // If the first slot is empty, place incoming there (new side).
          existing.first = incoming;
        } else {
          // Both slots are filled. Only upgrade missing owner info.
          // (Do not overwrite an existing owner assignment.)
          if (existing.first.index_unique < 0 && incoming.index_unique >= 0) {
            existing.first.index_unique = incoming.index_unique;
            existing.first.rank = incoming.rank;
          } else if (existing.second.index_unique < 0 && incoming.index_unique >= 0) {
            existing.second.index_unique = incoming.index_unique;
            existing.second.rank = incoming.rank;
          }
        }
      };
 
      // Merge both sides from the received pair.
      ingest(rbuff_value[i].first);
      ingest(rbuff_value[i].second);
    } else {
      // First time we see this key on the receiving rank.
      rmap.insert({rbuff_key[i], rbuff_value[i]});
    }
  }
 
  for (size_t i = 0; i < nrecv; i++ ) {
    auto it = rmap.find(rbuff_key[i]);
    if(it != rmap.end()) rbuff_value[i] = it->second;
  }
 
  // sending the assigned data back to original rank
  grph.transpose();
  MPI_Exchange(grph, rbuff_value, ibuff_value, comm);
 
  for (size_t i = 0; i < dsize; i++ ) {
    auto it = map.find(sbuff_key[i]);
    if (it != map.end()) {
      // Overwrite with merged entry from exchange (rmap already combined both sides)
      it->second = ibuff_value[i];
    } else {
      // Insert new entry received from another rank
      // This ensures non-owner ranks get the complete mapping data
      map.insert({sbuff_key[i], ibuff_value[i]});
    }
  }
}
 
template<class K> inline
void assign_counter_vertices_tuple(hashmap::unordered_map<K, intersection_data>& map, const MPI_Comm comm)
{
    size_t dsize = map.size();
    vector<K> key_vec(dsize);
    vector<intersection_data> value_vec(dsize);
 
    size_t idx = 0;
    for (const auto& v : map) {
        key_vec[idx] = v.first;
        value_vec[idx] = v.second;
        idx++;
    }
 
    vector<int> perm, dest;
    emi_select_merge_destinations(key_vec, dsize, sizeof(K) + sizeof(intersection_data), dest, comm,
                                  "assign_counter_vertices_tuple",
                                  [](const K& key) { return hashmap::hash_ops<K>::hash(key); });
 
    commgraph<size_t> grph;
    grph.configure(dest, comm);
    size_t nrecv = sum(grph.rcnt);
 
    interval(perm, 0, dsize);
    binary_sort_copy(dest, perm);
 
    // fill send buffer and communicate
    vector<K> sbuff_key(dsize), rbuff_key(nrecv);
    vector<intersection_data> sbuff_value(dsize), ibuff_value(dsize), rbuff_value(nrecv);
    for (size_t i = 0; i < dsize; i++) {
        sbuff_key[i] = key_vec[perm[i]];
        sbuff_value[i] = value_vec[perm[i]];
    }
    MPI_Exchange(grph, sbuff_key, rbuff_key, comm);
    MPI_Exchange(grph, sbuff_value, rbuff_value, comm);
 
    hashmap::unordered_map<K, intersection_data> rmap;
    for (size_t i = 0; i < nrecv; i++) {
        auto it = rmap.find(rbuff_key[i]);
        if (it != rmap.end()) {
            intersection_data& map_val = it->second;
            intersection_data& r_val = rbuff_value[i];
 
            int first_index = -1;
            for (int j = 0; j < MAX_INTERSECTIONS; ++j) {
                if (map_val.tags[j] == -1 && map_val.ranks[j] == -1) {
                    first_index = j;
                    break;
                }
            }
 
            if (first_index == -1) continue; // This means all possible intersections have already been assigned. 
 
            for (int j = 0; j < MAX_INTERSECTIONS; ++j) {
                if (r_val.tags[j] == -1) continue;
 
                int rtag = r_val.tags[j];
                int rdata = r_val.data[j];
                int rrank = r_val.ranks[j];
 
                bool exist_rtag = false;
                for (int k = 0; k < first_index; ++k) {
                    if (map_val.tags[k] == rtag && map_val.ranks[k] == rrank) {
                        exist_rtag = true;
                        break;
                    }
                }
 
                if (!exist_rtag && first_index < MAX_INTERSECTIONS) {
                    map_val.tags[first_index] = rtag;
                    map_val.data[first_index] = rdata;
                    map_val.ranks[first_index] = rrank;
                    first_index++;
                }
            }
        }
        else {
            rmap.insert({ rbuff_key[i], rbuff_value[i] });
        }
    }
 
    for (size_t i = 0; i < nrecv; i++) {
        auto it = rmap.find(rbuff_key[i]);
        if (it != rmap.end()) rbuff_value[i] = it->second;
    }
 
    // sending the assigned data back to original rank
    grph.transpose();
    MPI_Exchange(grph, rbuff_value, ibuff_value, comm);
 
    for (size_t i = 0; i < dsize; i++) {
        auto it = map.find(sbuff_key[i]);
        if (it != map.end()) it->second = ibuff_value[i];
    }
}
 
template<class K, class V> inline
void assign_dof_on_counter_face(hashmap::unordered_map<K, V>  & map, const MPI_Comm comm)
{
  size_t dsize = map.size();
  vector<K>    key_vec   (dsize);
  vector<V>    value_vec (dsize);
 
  // make key and value vector of elements without counterparts
  // which are from different ranks
  size_t idx=0;
  for(const auto & v : map) {
      key_vec[idx]   = v.first;
      value_vec[idx] = v.second;
      idx++;
  }
  dsize = idx;
  key_vec.resize(dsize);
  value_vec.resize(dsize);
  
  vector<int> perm, dest;
  emi_select_merge_destinations(key_vec, dsize, sizeof(K) + sizeof(V), dest, comm,
                                "assign_dof_on_counter_face",
                                [](const K& key) { return key.first; });
 
  commgraph<size_t> grph;
  grph.configure(dest, comm);
  size_t nrecv = sum(grph.rcnt);
 
  interval(perm, 0, dsize);
  binary_sort_copy(dest, perm);
 
  // fill send buffer and communicate
  vector<K> sbuff_key(dsize),   rbuff_key(nrecv);
  vector<V> sbuff_value(dsize), ibuff_value(dsize), rbuff_value(nrecv);
  for(size_t i=0; i<dsize; i++) {
    sbuff_key[i]   = key_vec[perm[i]];
    sbuff_value[i] = value_vec[perm[i]];
  }
  MPI_Exchange(grph, sbuff_key, rbuff_key, comm);
  MPI_Exchange(grph, sbuff_value, rbuff_value, comm);
 
  hashmap::unordered_map<K,V>  rmap;
 
  for (size_t i = 0; i < nrecv; i++) {
    auto it = rmap.find(rbuff_key[i]);
    if(it != rmap.end()) {
      if(rbuff_value[i] != -1 and it->second == -1){
        it->second    = rbuff_value[i];// add new indices
      }
    } else {
      rmap.insert({rbuff_key[i], rbuff_value[i]});
    }
  }
 
  for (size_t i = 0; i < nrecv; i++ ) {
    auto it = rmap.find(rbuff_key[i]);
    if(it != rmap.end()) {
      if(rbuff_value[i] == -1)
        rbuff_value[i] = it->second;
    }
  }
 
  // sending the assigned data back to original rank
  grph.transpose();
  MPI_Exchange(grph, rbuff_value, ibuff_value, comm);
 
  for (size_t i = 0; i < dsize; i++ ) {
    auto it = map.find(sbuff_key[i]);
    if (it != map.end()) it->second = ibuff_value[i];
  }
}
 
template<class K, class V> inline
void assign_petsc_on_counter_face(hashmap::unordered_map<K, V>  & map, const MPI_Comm comm)
{
  size_t dsize = map.size();
  vector<K>    key_vec   (dsize);
  vector<V>    value_vec (dsize);
 
  // make key and value vector of elements without counterparts
  // which are from different ranks
  size_t idx=0;
  for(const auto & v : map) {
      key_vec[idx]   = v.first;
      value_vec[idx] = v.second;
      idx++;
  }
  dsize = idx;
  key_vec.resize(dsize);
  value_vec.resize(dsize);
  
  vector<int> perm, dest;
  emi_select_merge_destinations(key_vec, dsize, sizeof(K) + sizeof(V), dest, comm,
                                "assign_petsc_on_counter_face",
                                [](const K& key) { return key.first; });
 
  commgraph<size_t> grph;
  grph.configure(dest, comm);
  size_t nrecv = sum(grph.rcnt);
 
  interval(perm, 0, dsize);
  binary_sort_copy(dest, perm);
 
  // fill send buffer and communicate
  vector<K> sbuff_key(dsize),   rbuff_key(nrecv);
  vector<V> sbuff_value(dsize), ibuff_value(dsize), rbuff_value(nrecv);
  for(size_t i=0; i<dsize; i++) {
    sbuff_key[i]   = key_vec[perm[i]];
    sbuff_value[i] = value_vec[perm[i]];
  }
  MPI_Exchange(grph, sbuff_key, rbuff_key, comm);
  MPI_Exchange(grph, sbuff_value, rbuff_value, comm);
 
  hashmap::unordered_map<K,V>  rmap;
 
  for (size_t i = 0; i < nrecv; i++) {
    auto it = rmap.find(rbuff_key[i]);
    if(it != rmap.end()) {
      if(rbuff_value[i].second != -1 and it->second.second == -1){
        it->second    = rbuff_value[i];// add new indices
      }
    } else {
      rmap.insert({rbuff_key[i], rbuff_value[i]});
    }
  }
 
  for (size_t i = 0; i < nrecv; i++ ) {
    auto it = rmap.find(rbuff_key[i]);
    if(it != rmap.end()) {
      if(rbuff_value[i].second == -1)
        rbuff_value[i].second = it->second.second;
    }
  }
 
  // sending the assigned data back to original rank
  grph.transpose();
  MPI_Exchange(grph, rbuff_value, ibuff_value, comm);
 
  for (size_t i = 0; i < dsize; i++ ) {
    auto it = map.find(sbuff_key[i]);
    if (it != map.end()) it->second = ibuff_value[i];
  }
}
 
template<class T>
void sort_surf_local_indices(tuple<T> & ref, tuple<T> & target)
{
  vector<T> buff(2);
 
  buff[0] = ref.v1, buff[1] = ref.v2;
  binary_sort(buff);
 
  target.v1 = buff[0], target.v2 = buff[1];
}
 
template<class T>
void sort_surf_local_indices(triple<T> & ref, triple<T> & target)
{
  vector<T> buff(3);
 
  buff[0] = ref.v1, buff[1] = ref.v2, buff[2] = ref.v3;
  binary_sort(buff);
 
  target.v1 = buff[0], target.v2 = buff[1], target.v3 = buff[2];
}
 
template<class T>
void sort_surf_local_indices(quadruple<T> & ref, quadruple<T> & target)
{
  vector<T> buff(4);
 
  buff[0] = ref.v1, buff[1] = ref.v2, buff[2] = ref.v3, buff[3] = ref.v4;
  binary_sort(buff);
 
  target.v1 = buff[0], target.v2 = buff[1], target.v3 = buff[2], target.v4 = buff[3];
}
 
template<class K, class V> inline
void insert_surf_based_Tag(V & ref,
                           hashmap::unordered_map<K, std::pair<V, V>> & surfmap)
{
  K surf;
 
  sort_surf_local_indices(ref.points, surf);
 
  auto it = surfmap.find(surf);
  if (it != surfmap.end()) {
    // Add as a counter face only if this face exists in surfMap and its tag number
    // differs from the tag number of the first saved pair.
    if ( it->second.first.tag != ref.tag )
      it->second.second = ref;
  } else {
    std::pair<V, V> face;
    face.first = ref;
    surfmap.insert({surf,face});
  }
}
 
template<class T, class W, class V, class U> inline
void insert_surf_emi( int rank, 
                      SF::vector<T> const & ref_con, 
                      SF::vector<T> const & ptsData,
                      const T eidx, const T tag,
                      const std::vector<T> & line_con, const std::vector<T> & surf_con, const std::vector<T> & qsurf_con,
                      hashmap::unordered_map<tuple<T>, std::pair<W, W> > & line_face,
                      hashmap::unordered_map<triple<T>, std::pair<V, V> > & surfmap,
                      hashmap::unordered_map<quadruple<T>, std::pair<U, U> > & qsurfmap)
{
  W line;
  line.eidx = eidx;
  line.tag  = tag;
  V face;
  face.eidx = eidx;
  face.tag  = tag;
  U qface;
  qface.eidx = eidx;
  qface.tag  = tag;  
 
  // Mark faces as potential candidates if all of their vertices lie on the interface.
  for (size_t i = 0; i < line_con.size(); i += 2) {
    if (ptsData[line_con[i    ] - 1] > 0 &&
        ptsData[line_con[i + 1] - 1] > 0 ) {
      line.points.v1 = ref_con[line_con[i    ] - 1]; 
      line.points.v2 = ref_con[line_con[i + 1] - 1]; 
      line.rank = rank;
      insert_surf_based_Tag(line, line_face);
    }
  }
 
  for (size_t i = 0; i < surf_con.size(); i += 3) {
    if (ptsData[surf_con[i    ] - 1] > 0 &&
        ptsData[surf_con[i + 1] - 1] > 0 &&
        ptsData[surf_con[i + 2] - 1] > 0 ) {
      face.points.v1 = ref_con[surf_con[i    ] - 1]; 
      face.points.v2 = ref_con[surf_con[i + 1] - 1]; 
      face.points.v3 = ref_con[surf_con[i + 2] - 1];
      face.rank = rank; 
      insert_surf_based_Tag(face, surfmap);
    }
  }
 
  for (size_t i = 0; i < qsurf_con.size(); i += 4) {
    if (ptsData[qsurf_con[i    ] - 1] > 0 &&
        ptsData[qsurf_con[i + 1] - 1] > 0 &&
        ptsData[qsurf_con[i + 2] - 1] > 0 &&
        ptsData[qsurf_con[i + 3] - 1] > 0) {
      qface.points.v1 = ref_con[qsurf_con[i    ] - 1]; 
      qface.points.v2 = ref_con[qsurf_con[i + 1] - 1]; 
      qface.points.v3 = ref_con[qsurf_con[i + 2] - 1]; 
      qface.points.v4 = ref_con[qsurf_con[i + 3] - 1]; 
      qface.rank = rank; 
      insert_surf_based_Tag(qface, qsurfmap);
    }
  }
}
 
template<class T, class V>
struct emi_face {
  V points;
  T eidx = -1;
  T tag  = 0;
  T rank = -1;
  T mem = 1; // default the face on membrane, otherwise will be asssign to 2 for gap-junction face
  bool mark_to_take = false;  // if the face is selected on one of the ranks on the surface mesh with unique face, it will be true
                              // we assume that the face will be selected on the smaller rank. If both has the same rank, we choose the first pair.
  T index_unique = -1;        // index of face on emi_surface_unique_face_msh
  T index_one = -1;           // index of face on emi_surface_msh
};
 
template<class T>
struct emi_unique_face {
  T index_unique = -1;
  T index_one = -1;
  T rank = -1;
};
 
template<class T>
struct emi_index_rank {
  T index = -1;
  T rank = -1;
};
template<class T, class S> inline
void compute_surface_with_tags(meshdata<T,S> & mesh, 
                               const SF_nbr numbering,
                               hashmap::unordered_map<T,T> & vertex2ptsdata,
                               hashmap::unordered_set<int> & extra_tags,
                               hashmap::unordered_map<tuple<T>, 
                               std::pair<emi_face<T,tuple<T>>, emi_face<T,tuple<T>>>> & line_face,
                               hashmap::unordered_map<triple<T>, 
                               std::pair<emi_face<T,triple<T>>, emi_face<T,triple<T>>>>     & tri_surf,
                               hashmap::unordered_map<quadruple<T>, 
                               std::pair<emi_face<T,quadruple<T>>, emi_face<T,quadruple<T>>>> & quad_surf)
{
  MPI_Comm comm = SF_COMM;
  int      size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  const T* con = mesh.con.data();
  const T* nbr = mesh.get_numbering(numbering).data();
  const SF::vector<mesh_int_t> & rnod = mesh.get_numbering(numbering);
 
  hashmap::unordered_map<mesh_int_t,mesh_int_t> g2ptsData;
 
  for(size_t i=0; i<mesh.con.size(); i++){
      g2ptsData[rnod[con[i]]] = vertex2ptsdata[rnod[con[i]]];
  }
 
  const vector<T> & ref_eidx = mesh.get_numbering(NBR_ELEM_REF);
 
  // Collect potential faces whose three vertices lie on the EMI interfaces.
  // The ptsData, computed earlier from the input mesh, is used to check
  // the intersection of each vertex with different tag numbers.
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    T tag = mesh.tag[eidx];
    T size_elem = mesh.dsp[eidx+1]-mesh.dsp[eidx];
    
    vector<T> dofvec(size_elem); // reference nodes of each element
    vector<T> ptsDatavec(size_elem); // point data on the nodes of each element
 
    for (int n = mesh.dsp[eidx], i = 0; n < mesh.dsp[eidx+1];n++,i++)
    {
      dofvec[i] = rnod[con[n]];
      ptsDatavec[i] =  g2ptsData[rnod[con[n]]];
    }
    std::vector<T> surf_con ;
    std::vector<T> qsurf_con;
    std::vector<T> line_con;
    switch(mesh.type[eidx]) {
      // I need to find the right order of lines in 2D 
 
      // I guess that in the 2D case we would define the interfaces as
      // lines instead of surfaces.
    case Tri: {
      line_con = {1, 2,
                  2, 3,
                  3, 1};
      surf_con  = {};
      qsurf_con  = {};
      break;
    }
 
    case Quad: {
      line_con = {1, 2,
                  2, 3,
                  3, 4,
                  4, 1};
      surf_con  = {};
      qsurf_con  = {};
      break;
    }
 
    case Tetra: {
      // surfaces are (2,3,1) , (1,4,2) , (2,4,3) , (1,3,4)
      line_con = {};
      surf_con = {2,3,1,
                  1,4,2, 
                  2,4,3,
                  1,3,4};
      qsurf_con = {};
      break; 
    }
    case Pyramid: {
      // surfaces are (1,5,2) , (2,5,3) , (3,5,4) , (4,5,1) , (1,2,3,4)
      line_con = {};
      surf_con = {1,5,2,
                  2,5,3, 
                  3,5,4,
                  4,5,1};
      qsurf_con = {1,2,3,4};
      break;
    }
    case Prism: {
      // surfaces are (1,2,3) , (4,5,6) , (1,2,6,4) , (2,3,5,6) , (3,1,4,5)
      line_con = {};
      surf_con = {1,2,3, 
                  4,5,6};
      qsurf_con = {1,2,6,4,
                   2,3,5,6, 
                   3,1,4,5};
      break;
    }
    case Hexa: {
      // surfaces are (1,2,3,4) , (3,2,8,7) , (4,3,7,6) , (1,4,6,5) , (2,1,5,8), (6,7,8,5)
      line_con = {};
      surf_con  = {};
      qsurf_con = {1,2,3,4,
                   3,2,8,7,
                   4,3,7,6,
                   1,4,6,5,
                   2,1,5,8,
                   6,7,8,5 };
      break;
    }
    default:
      fprintf(stderr, "%s error: Unsupported element in surface computation!\n", __func__);
      exit(1);
    }
    insert_surf_emi(rank, dofvec, ptsDatavec, ref_eidx[eidx], mesh.tag[eidx], line_con, surf_con, qsurf_con, line_face, tri_surf, quad_surf);
  }
 
  // Find the counter faces corresponding to the selected potential faces.
  assign_counter_face(line_face, mesh.comm);
  assign_counter_face(tri_surf, mesh.comm);
  assign_counter_face(quad_surf, mesh.comm);
 
  // Removed potential faces without matching counter faces and 
  // Removed potential faces where both associated tags belong to the extracellular domain. 
  for(auto it = line_face.begin(); it != line_face.end(); ) {
    if( it->second.second.eidx == -1 ||(extra_tags.find(it->second.first.tag) != extra_tags.end() &&  extra_tags.find(it->second.second.tag) != extra_tags.end()))
      it = line_face.erase(it);
    else
      ++it;
  }
 
  for(auto it = tri_surf.begin(); it != tri_surf.end(); ) {
    if( it->second.second.eidx == -1 ||(extra_tags.find(it->second.first.tag) != extra_tags.end() &&  extra_tags.find(it->second.second.tag) != extra_tags.end()))
      it = tri_surf.erase(it);
    else
      ++it;
  }
 
  for(auto it = quad_surf.begin(); it != quad_surf.end(); ) {
    if( it->second.second.eidx == -1 ||(extra_tags.find(it->second.first.tag) != extra_tags.end() &&  extra_tags.find(it->second.second.tag) != extra_tags.end()))
      it = quad_surf.erase(it);
    else
      ++it;
  }
}
 
inline bool should_take_first(int tag1, int tag2, 
                              const hashmap::unordered_set<int>& extra_tags,
                              const hashmap::unordered_set<int>& intra_tags)
{
  bool tag1_is_intra = (intra_tags.find(tag1) != intra_tags.end());
  bool tag2_is_intra = (intra_tags.find(tag2) != intra_tags.end());
  bool tag1_is_extra = (extra_tags.find(tag1) != extra_tags.end());
  bool tag2_is_extra = (extra_tags.find(tag2) != extra_tags.end());
  
  // Case 1: Both are intra tags (gap junction) - take smaller tag
  if(tag1_is_intra && tag2_is_intra) {
    return tag1 < tag2;
  }
  
  // Case 2: One is extra, one is intra (membrane) - take intra
  if(tag1_is_intra && tag2_is_extra) {
    return true;  // take tag1 (intra)
  }
  if(tag2_is_intra && tag1_is_extra) {
    return false; // take tag2 (intra)
  }
  
  // Default fallback: take smaller tag
  return tag1 < tag2;
}
 
inline
void assign_map_between_elem_oneface_and_elem_uniqueFace(int rank, 
                                                emi_unique_face <mesh_int_t> first,
                                                emi_unique_face <mesh_int_t> second,
                                                hashmap::unordered_map<mesh_int_t, std::pair<emi_index_rank<mesh_int_t>, emi_index_rank<mesh_int_t>>> & map_elem_uniqueFace_to_elem_oneface,
                                                hashmap::unordered_map<mesh_int_t, emi_index_rank<mesh_int_t>> & map_elem_oneface_to_elem_uniqueFace)
{
  // Case 1: Both face and counter face belong to the same rank
  if(second.index_one==-1 && second.index_unique==-1){
    second.index_one = first.index_one;
    second.index_unique = first.index_unique;
  }
  // Case 2: Cross-rank faces - set up mapping for non-owner ranks
  else if(first.rank != second.rank) {
    // Find our local index_one (from whichever field matches our rank)
    int our_index_one = (first.rank == rank) ? first.index_one : 
                        (second.rank == rank) ? second.index_one : -1;
    // Find owner's index_unique (from whichever field has valid index_unique)
    int owner_index_unique = (first.index_unique >= 0) ? first.index_unique :
                             (second.index_unique >= 0) ? second.index_unique : -1;
    int owner_rank = (first.index_unique >= 0) ? first.rank : 
                     (second.index_unique >= 0) ? second.rank : -1;
    
    // Only set mapping for non-owner ranks (owner already set it in STEP 3)
    if(our_index_one >= 0 && owner_index_unique >= 0 && owner_rank >= 0 && owner_rank != rank) {
      map_elem_oneface_to_elem_uniqueFace[our_index_one].index = owner_index_unique;
      map_elem_oneface_to_elem_uniqueFace[our_index_one].rank = owner_rank;
    }
  }
 
  // Populate map_elem_uniqueFace_to_elem_oneface (only on owner rank where index_unique >= 0)
  // Checked both first and second to find which one represents this rank as owner
  // After MPI exchange, the owner info could be in either field
  int local_index_unique = -1;
  emi_index_rank<mesh_int_t> value1, value2;
  
  if (first.index_unique >= 0 && first.rank == rank) {
    // This rank owns the unique face, info is in 'first'
    local_index_unique = first.index_unique;
    value1.index = first.index_one;
    value1.rank = first.rank;
    value2.index = second.index_one;
    value2.rank = second.rank;
  } else if (second.index_unique >= 0 && second.rank == rank) {
    // This rank owns the unique face, info is in 'second'
    local_index_unique = second.index_unique;
    value1.index = second.index_one;
    value1.rank = second.rank;
    value2.index = first.index_one;
    value2.rank = first.rank;
  }
  
  if (local_index_unique >= 0) {
    auto itr = map_elem_uniqueFace_to_elem_oneface.find(local_index_unique);
    if (itr == map_elem_uniqueFace_to_elem_oneface.end()) {
      std::pair <emi_index_rank<mesh_int_t>,emi_index_rank<mesh_int_t>> value;
      value = std::make_pair(value1, value2);
      map_elem_uniqueFace_to_elem_oneface.insert({local_index_unique, value});
    }
  }
}
 
template <class T, class Key, class PointsKey>
inline void assign_ownership_rank_on_faces( const Key& key,
                                            std::pair<emi_face<T, PointsKey>, emi_face<T, PointsKey>>& v,  // (face, counter-face)
                                            int rank,
                                            mesh_int_t idx_oneface,
                                            mesh_int_t& lsize_unique,
                                            const hashmap::unordered_set<int>& extra_tags,
                                            const hashmap::unordered_set<int>& intra_tags,
                                            hashmap::unordered_map<Key,std::pair<emi_unique_face<mesh_int_t>, emi_unique_face<mesh_int_t>>>& unique_face_to_elements,
                                            hashmap::unordered_map<mesh_int_t, emi_index_rank<mesh_int_t>>& map_elem_oneface_to_elem_uniqueFace)
{
  const bool same_rank = (v.first.rank == v.second.rank);
 
  // Decide owner rank in an order-independent way:
  // - If one side is intra, always take the intra side
  // - If both intra, take smaller tag (tie -> smaller rank)
  // - Otherwise, take smaller tag (tie -> smaller rank)
  const bool first_intra  = (intra_tags.find(v.first.tag)  != intra_tags.end());
  const bool second_intra = (intra_tags.find(v.second.tag) != intra_tags.end());
 
  // Decide which side to take based purely on tags/intra (order-independent)
  bool take_first_by_tags = true;
  if (first_intra != second_intra) {
    take_first_by_tags = first_intra; // take intra side
  } else if (v.first.tag != v.second.tag) {
    take_first_by_tags = (v.first.tag < v.second.tag);
  } else {
    // Same tag (should be rare for interfaces); fall back to smaller rank to be deterministic
    take_first_by_tags = (v.first.rank <= v.second.rank);
  }
 
  int owner_rank = take_first_by_tags ? v.first.rank : v.second.rank;
 
  const bool i_am_owner = (rank == owner_rank);
 
  // If I'm the owner, mark the chosen side and set mapping now.
  if (i_am_owner) {
    auto& chosen = take_first_by_tags ? v.first : v.second;
 
    chosen.mark_to_take = true;
    chosen.index_unique = lsize_unique;
    chosen.index_one    = idx_oneface;
 
    map_elem_oneface_to_elem_uniqueFace[idx_oneface].index = lsize_unique;
    map_elem_oneface_to_elem_uniqueFace[idx_oneface].rank  = rank;
  }
 
  // Insert per-face record once:
  // first = "my side if I'm owner else -1" + rank info
  // second = placeholder for the other rank (or same rank if same_rank)
  auto it = unique_face_to_elements.find(key);
  if (it == unique_face_to_elements.end()) {
    emi_unique_face<mesh_int_t> first;
    emi_unique_face<mesh_int_t> second;
 
    if (i_am_owner) {
      // Owner knows its local unique index.
      first.index_unique = lsize_unique;
      first.index_one    = idx_oneface;
      first.rank         = rank;
 
      // "other side" rank:
      second.index_unique = -1;
      second.index_one    = -1;
      second.rank         = same_rank ? rank : (rank == v.first.rank ? v.second.rank : v.first.rank);
    } else {
      // Non-owner does not know unique index yet; MPI step will fill it.
      first.index_unique = -1;
      first.index_one    = idx_oneface;
      first.rank = rank;
 
      second.index_unique = -1;
      second.index_one    = -1;
      second.rank         = owner_rank; // this is the rank that will eventually provide index_unique
    }
 
    unique_face_to_elements.insert({ key, std::make_pair(first, second) });
  } else {
    // Update existing entry to avoid losing owner info when non-owner inserted first.
    auto& existing = it->second;
 
    // Choose a slot to update for this rank.
    emi_unique_face<mesh_int_t>* slot = nullptr;
    if (existing.first.rank == rank || existing.first.rank < 0) {
      slot = &existing.first;
    } else if (existing.second.rank == rank || existing.second.rank < 0) {
      slot = &existing.second;
    } else if (existing.first.index_one == idx_oneface) {
      slot = &existing.first;
    } else if (existing.second.index_one == idx_oneface) {
      slot = &existing.second;
    } else if (existing.first.index_unique < 0) {
      slot = &existing.first;
    } else if (existing.second.index_unique < 0) {
      slot = &existing.second;
    }
 
    if (slot) {
      if (slot->rank < 0) slot->rank = rank;
      if (slot->index_one < 0) slot->index_one = idx_oneface;
      if (i_am_owner && slot->index_unique < 0) {
        slot->index_unique = lsize_unique;
      }
    }
  }
 
  // Only the owning rank increments lsize_unique
  if (i_am_owner) {
    lsize_unique += 1;
  }
}
 
 
template<class T, class S> inline
void extract_face_based_tags(meshdata<T,S> & mesh, 
                                     const SF_nbr numbering,
                                     hashmap::unordered_map<T,T> & vertex2ptsdata, 
                                     hashmap::unordered_map<tuple<T> , 
                                     std::pair<emi_face<T,tuple<T>>, emi_face<T,tuple<T>>>> & line_face,
                                     hashmap::unordered_map<triple<T>, 
                                     std::pair<emi_face<T,triple<T>>, emi_face<T,triple<T>>>> & tri_face,
                                     hashmap::unordered_map<quadruple<T>, 
                                     std::pair<emi_face<T,quadruple<T>>, emi_face<T,quadruple<T>>>> & quad_face,
                                     hashmap::unordered_set<int> & extra_tags,
                                     hashmap::unordered_set<int> & intra_tags,
                                     meshdata<T,S> & surfmesh,
                                     meshdata<T,S> & surfmesh_w_counter,
                                     meshdata<T,S> & surfmesh_unique_face,
                                     hashmap::unordered_map<mesh_int_t, std::pair<emi_index_rank<mesh_int_t>, emi_index_rank<mesh_int_t>>> & map_elem_uniqueFace_to_elem_oneface,
                                     hashmap::unordered_map<mesh_int_t, emi_index_rank<mesh_int_t>> & map_elem_oneface_to_elem_uniqueFace)
{
  // ============================================================================
  // STEP 1: Initialize surface meshes and extract interface faces
  // ============================================================================
  surfmesh.register_numbering(SF::NBR_REF);
  surfmesh_w_counter.register_numbering(SF::NBR_REF);
  surfmesh_unique_face.register_numbering(SF::NBR_REF); 
 
  // Extract faces with counter faces on the interfaces to build surface mesh
 
  compute_surface_with_tags(mesh, numbering,  vertex2ptsdata, extra_tags, line_face, tri_face, quad_face);
 
  int      size, rank;
  MPI_Comm_size(mesh.comm, &size);
  MPI_Comm_rank(mesh.comm, &rank);
 
  long int g_num_line = line_face.size();
  long int g_num_tri  = tri_face.size();
  long int g_num_quad = quad_face.size();
  MPI_Allreduce(MPI_IN_PLACE, &g_num_line,  1, MPI_LONG, MPI_SUM, mesh.comm);
  MPI_Allreduce(MPI_IN_PLACE, &g_num_tri,  1, MPI_LONG, MPI_SUM, mesh.comm);
  MPI_Allreduce(MPI_IN_PLACE, &g_num_quad, 1, MPI_LONG, MPI_SUM, mesh.comm);
 
  // Warn if any extra-extra faces remain (should be none after filtering)
  {
    long int local_extra_extra = 0;
    for (const auto& [key, v] : line_face) {
      if (extra_tags.find(v.first.tag) != extra_tags.end() &&
          extra_tags.find(v.second.tag) != extra_tags.end()) {
        local_extra_extra++;
      }
    }
    for (const auto& [key, v] : tri_face) {
      if (extra_tags.find(v.first.tag) != extra_tags.end() &&
          extra_tags.find(v.second.tag) != extra_tags.end()) {
        local_extra_extra++;
      }
    }
    for (const auto& [key, v] : quad_face) {
      if (extra_tags.find(v.first.tag) != extra_tags.end() &&
          extra_tags.find(v.second.tag) != extra_tags.end()) {
        local_extra_extra++;
      }
    }
    long int global_extra_extra = 0;
    MPI_Allreduce(&local_extra_extra, &global_extra_extra, 1, MPI_LONG, MPI_SUM, mesh.comm);
    if (global_extra_extra > 0 && rank == 0) {
      fprintf(stderr, "WARN: extra-extra faces detected before filtering: %ld\n", global_extra_extra);
    }
  }
 
  surfmesh.g_numelem = g_num_line + g_num_tri + g_num_quad;
  surfmesh.l_numelem = line_face.size() + tri_face.size() + quad_face.size();
 
  vector<T> cnt(surfmesh.l_numelem, 0);
  surfmesh.tag.resize(surfmesh.l_numelem);
  surfmesh.type.resize(surfmesh.l_numelem);
  surfmesh.con.resize(line_face.size() * 2 + tri_face.size() * 3 + quad_face.size() * 4);
 
  // ============================================================================
  // STEP 2: Sort face keys for deterministic ordering across MPI ranks
  // ============================================================================
  // To ensure deterministic order, we must sort the keys of the hashmaps before iterating.
  // This is crucial for MPI consistency - all ranks must process faces in the same order.
  vector<tuple<T>> line_keys;
  for(auto const& [key, val] : line_face) line_keys.push_back(key);
  std::sort(line_keys.begin(), line_keys.end());
 
  vector<triple<T>> tri_keys;
  for(auto const& [key, val] : tri_face) tri_keys.push_back(key);
  std::sort(tri_keys.begin(), tri_keys.end());
 
  vector<quadruple<T>> quad_keys;
  for(auto const& [key, val] : quad_face) quad_keys.push_back(key);
  std::sort(quad_keys.begin(), quad_keys.end());
  
  // Temporary hashmaps to track unique-face assignments before MPI communication
  hashmap::unordered_map<SF::tuple<mesh_int_t>, std::pair<emi_unique_face <mesh_int_t>,emi_unique_face<mesh_int_t>>> line_unique_face_to_elements;
  hashmap::unordered_map<SF::triple<mesh_int_t>, std::pair<emi_unique_face<mesh_int_t>,emi_unique_face<mesh_int_t>>> tri_unique_face_to_elements;
  hashmap::unordered_map<SF::quadruple<mesh_int_t>, std::pair<emi_unique_face<mesh_int_t>,emi_unique_face<mesh_int_t>>> quad_unique_face_to_elements;
 
  // Counters for local element indices:
  // - idx: index into surfmesh (one-side mesh)
  // - lsize_*: count of elements in surfmesh_w_counter
  // - lsize_unique_*: count of elements in surfmesh_unique_face
  mesh_int_t idx = 0, cidx = 0;
  mesh_int_t lsize_line = 0;
  mesh_int_t lsize_tri = 0;
  mesh_int_t lsize_quad = 0;
 
  mesh_int_t lsize_unique_line = 0;
  mesh_int_t lsize_unique_tri = 0;
  mesh_int_t lsize_unique_quad = 0;
 
  // ============================================================================
  // STEP 3: Process each face and determine unique-face ownership
  // ============================================================================
  // For each interface face, we:
  // 1. Add it to surfmesh
  for(auto const& key : line_keys) {
    auto& v = line_face.at(key);
    const bool local_involved = (v.first.rank == rank) || (v.second.rank == rank);
    if (!local_involved) {
      continue;
    }
    cnt[idx] = 2;
 
    // Mark gap-junction faces explicitly; membrane is the default.
    if ((intra_tags.find(v.first.tag) != intra_tags.end()) && (intra_tags.find(v.second.tag) != intra_tags.end())) {
      v.first.mem = 2; // to set the face located on gap-junction
      v.second.mem = 2; // to set the face located on gap-junction
    }
 
    // changing the strategy to select the face with the one with the same rank 
    emi_face<T,tuple<T>> surf_neighbor = 
      (v.first.rank == rank) ? v.first : v.second;
 
    surfmesh.type[idx]     = Line;
    surfmesh.tag[idx]      = surf_neighbor.tag;
    surfmesh.con[cidx + 0] = surf_neighbor.points.v1;
    surfmesh.con[cidx + 1] = surf_neighbor.points.v2;
 
    // If the ranks are equal, both faces will be added later.
    if(v.first.rank == v.second.rank)
      lsize_line+=2;
    else
      lsize_line+=1;
   
    assign_ownership_rank_on_faces<T, tuple<T>, tuple<T>>( key, v, rank,
                                                           idx,
                                                           lsize_unique_line,
                                                           extra_tags, intra_tags,
                                                           line_unique_face_to_elements,
                                                           map_elem_oneface_to_elem_uniqueFace);
    idx  += 1;
    cidx += 2;
  }
 
  for(auto const& key : tri_keys) {
    auto& v = tri_face.at(key);
    const bool local_involved = (v.first.rank == rank) || (v.second.rank == rank);
    if (!local_involved) {
      continue;
    }
    cnt[idx] = 3;
 
    if ((intra_tags.find(v.first.tag) != intra_tags.end()) && (intra_tags.find(v.second.tag) != intra_tags.end())) {
      v.first.mem = 2; // to set the face located on gap-junction
      v.second.mem = 2; // to set the face located on gap-junction
    }
 
    // changing the strategy to select the face with the one with the same rank 
    emi_face<T,triple<T>> surf_neighbor = 
      (v.first.rank == rank) ? v.first : v.second;
 
    surfmesh.type[idx]     = Tri;
    surfmesh.tag[idx]      = surf_neighbor.tag; 
    surfmesh.con[cidx + 0] = surf_neighbor.points.v1;
    surfmesh.con[cidx + 1] = surf_neighbor.points.v2;
    surfmesh.con[cidx + 2] = surf_neighbor.points.v3;
 
    // If the ranks are equal, both faces will be added later.
    if(v.first.rank == v.second.rank)
      lsize_tri+=2;
    else
      lsize_tri+=1;
 
    assign_ownership_rank_on_faces<T, triple<T>, triple<T>>(  key, v, rank,
                                                              idx,
                                                              lsize_unique_tri,
                                                              extra_tags, intra_tags,
                                                              tri_unique_face_to_elements,
                                                              map_elem_oneface_to_elem_uniqueFace);
    idx  += 1;
    cidx += 3;
  }
 
  for(auto const& key : quad_keys) {
    auto& v = quad_face.at(key);
    const bool local_involved = (v.first.rank == rank) || (v.second.rank == rank);
    if (!local_involved) {
      continue;
    }
    cnt[idx] = 4;
 
    if ((intra_tags.find(v.first.tag) != intra_tags.end()) && (intra_tags.find(v.second.tag) != intra_tags.end())) {
      v.first.mem = 2; // to set the face located on gap-junction
      v.second.mem = 2; // to set the face located on gap-junction
    }
    
    // changing the strategy to select the face with the one with the same rank
    emi_face<T,quadruple<T>> qsurf_neighbor = 
      (v.first.rank == rank) ? v.first : v.second;
 
    surfmesh.type[idx]     = Quad;
    surfmesh.tag[idx]      = qsurf_neighbor.tag; 
    surfmesh.con[cidx + 0] = qsurf_neighbor.points.v1;
    surfmesh.con[cidx + 1] = qsurf_neighbor.points.v2;
    surfmesh.con[cidx + 2] = qsurf_neighbor.points.v3;
    surfmesh.con[cidx + 3] = qsurf_neighbor.points.v4;
 
    // If the ranks are equal, both faces will be added later.
    if(v.first.rank == v.second.rank)
      lsize_quad+=2;
    else
      lsize_quad+=1;
 
    assign_ownership_rank_on_faces<T, quadruple<T>, quadruple<T>>(  key, v, rank,
                                                                    idx,
                                                                    lsize_unique_quad,
                                                                    extra_tags, intra_tags,
                                                                    quad_unique_face_to_elements,
                                                                    map_elem_oneface_to_elem_uniqueFace);
    idx  += 1;
    cidx += 4;
  }
 
  surfmesh.l_numelem = idx;
  surfmesh.tag.resize(idx);
  surfmesh.type.resize(idx);
  cnt.resize(idx);
  dsp_from_cnt(cnt, surfmesh.dsp);
  surfmesh.con.resize(cidx);
 
  long int g_num_surf = surfmesh.l_numelem;
  MPI_Allreduce(MPI_IN_PLACE, &g_num_surf, 1, MPI_LONG, MPI_SUM, mesh.comm);
  surfmesh.g_numelem = g_num_surf;
 
  if(rank ==0){
    std::cout << "surfmesh.g_numelem: " << surfmesh.g_numelem <<std::endl;
    std::cout << "surfmesh.l_numelem: " << surfmesh.l_numelem <<std::endl;
  }
 
  {
    long int g_num_w_counter_line = lsize_line;
    long int g_num_w_counter_tri  = lsize_tri;
    long int g_num_w_counter_quad = lsize_quad;
 
    MPI_Allreduce(MPI_IN_PLACE, &g_num_w_counter_line,  1, MPI_LONG, MPI_SUM, SF_COMM);
    MPI_Allreduce(MPI_IN_PLACE, &g_num_w_counter_tri,  1, MPI_LONG, MPI_SUM, SF_COMM);
    MPI_Allreduce(MPI_IN_PLACE, &g_num_w_counter_quad,  1, MPI_LONG, MPI_SUM, SF_COMM);
 
    // Assign global and local elements for the surface mesh used in ionic computation.
    surfmesh_w_counter.g_numelem = g_num_w_counter_line+g_num_w_counter_tri+g_num_w_counter_quad;
    surfmesh_w_counter.l_numelem = lsize_line + lsize_tri + lsize_quad;
 
    surfmesh_w_counter.tag.resize(surfmesh_w_counter.l_numelem);
    surfmesh_w_counter.type.resize(surfmesh_w_counter.l_numelem);
    surfmesh_w_counter.con.resize(lsize_line * 2 + lsize_tri * 3 + lsize_quad * 4);
 
    if(rank ==0){
      std::cout << "surfmesh_w_counter.g_numelem: " << surfmesh_w_counter.g_numelem <<std::endl;
      std::cout << "surfmesh_w_counter.l_numelem: " << surfmesh_w_counter.l_numelem <<std::endl;
    }
  }
 
  {
    long int g_num_w_unique_line = lsize_unique_line;
    long int g_num_w_unique_tri  = lsize_unique_tri;
    long int g_num_w_unique_quad = lsize_unique_quad;
 
    MPI_Allreduce(MPI_IN_PLACE, &g_num_w_unique_line,  1, MPI_LONG, MPI_SUM, SF_COMM);
    MPI_Allreduce(MPI_IN_PLACE, &g_num_w_unique_tri,  1, MPI_LONG, MPI_SUM, SF_COMM);
    MPI_Allreduce(MPI_IN_PLACE, &g_num_w_unique_quad,  1, MPI_LONG, MPI_SUM, SF_COMM);
 
    // Assign global and local elements for the surface mesh used in ionic computation.
    surfmesh_unique_face.g_numelem = g_num_w_unique_line+g_num_w_unique_tri+g_num_w_unique_quad;
    surfmesh_unique_face.l_numelem = lsize_unique_line + lsize_unique_tri + lsize_unique_quad;
 
    surfmesh_unique_face.tag.resize(surfmesh_unique_face.l_numelem);
    surfmesh_unique_face.type.resize(surfmesh_unique_face.l_numelem);
    surfmesh_unique_face.con.resize(lsize_unique_line * 2 + lsize_unique_tri * 3 + lsize_unique_quad * 4);
 
    if(rank ==0){
      std::cout << "surfmesh_unique_face.g_numelem: " << surfmesh_unique_face.g_numelem <<std::endl;
      std::cout << "surfmesh_unique_face.l_numelem: " << surfmesh_unique_face.l_numelem <<std::endl;
    }
  }
 
  // ============================================================================
  // STEP 4: MPI communication to share unique-face ownership between ranks
  // ============================================================================
  // When faces span multiple ranks, only one rank "owns" the unique-face.
  // assign_unique_first_face communicates the owner's index_unique to non-owner ranks
  // so they can set up their map_elem_oneface_to_elem_uniqueFace mappings correctly.
  
  // DEBUG: Before exchange
  #ifdef EMI_DEBUG_MESH
  {
    fprintf(stderr, "RANK %d BEFORE exchange: line_hash=%zu, tri_hash=%zu, quad_hash=%zu (expected unique total=%ld)\n",
          rank, line_unique_face_to_elements.size(), tri_unique_face_to_elements.size(), 
          quad_unique_face_to_elements.size(), lsize_unique_line + lsize_unique_tri + lsize_unique_quad);
    fflush(stderr);
  } 
  #endif
 
  if(size>1)
  {
    // Find the counter faces corresponding to the selected potential faces.
    assign_unique_first_face(line_unique_face_to_elements, mesh.comm);
    assign_unique_first_face(tri_unique_face_to_elements, mesh.comm);
    assign_unique_first_face(quad_unique_face_to_elements, mesh.comm);    
  }  
 
  // DEBUG: After exchange
  #ifdef EMI_DEBUG_MESH
  {
    int valid_line = 0, valid_tri = 0, valid_quad = 0;
    for (const auto& [key, val] : line_unique_face_to_elements) {
      if (val.first.index_unique >= 0 || val.second.index_unique >= 0) valid_line++;
    }
    for (const auto& [key, val] : tri_unique_face_to_elements) {
      if (val.first.index_unique >= 0 || val.second.index_unique >= 0) valid_tri++;
    }
    for (const auto& [key, val] : quad_unique_face_to_elements) {
      if (val.first.index_unique >= 0 || val.second.index_unique >= 0) valid_quad++;
    }
    fprintf(stderr, "RANK %d AFTER exchange: line_hash=%zu, tri_hash=%zu, quad_hash=%zu\n",
            rank, line_unique_face_to_elements.size(), tri_unique_face_to_elements.size(), 
            quad_unique_face_to_elements.size());
    fprintf(stderr, "RANK %d VALID entries (index_unique>=0): line=%d, tri=%d, quad=%d\n",
            rank, valid_line, valid_tri, valid_quad);
    fflush(stderr);
  }
  #endif
 
  // ============================================================================
  // STEP 5: Finalize mappings after MPI communication
  // ============================================================================
  for(auto it = line_unique_face_to_elements.begin(); it != line_unique_face_to_elements.end(); ++it) {
    
    auto& first = it->second.first;
    auto& second = it->second.second;
    
    assign_map_between_elem_oneface_and_elem_uniqueFace(rank, first, second, map_elem_uniqueFace_to_elem_oneface, map_elem_oneface_to_elem_uniqueFace);
  }
 
  for(auto it = tri_unique_face_to_elements.begin(); it != tri_unique_face_to_elements.end(); ++it) {
    auto& first = it->second.first;
    auto& second = it->second.second;
    
    assign_map_between_elem_oneface_and_elem_uniqueFace(rank, first, second, map_elem_uniqueFace_to_elem_oneface, map_elem_oneface_to_elem_uniqueFace);
  }
 
  for(auto it = quad_unique_face_to_elements.begin(); it != quad_unique_face_to_elements.end(); ++it) {
    auto& first = it->second.first;
    auto& second = it->second.second;
    
    assign_map_between_elem_oneface_and_elem_uniqueFace(rank, first, second, map_elem_uniqueFace_to_elem_oneface, map_elem_oneface_to_elem_uniqueFace);
  }
}
 
template<class T> inline
void complete_map_vertex_to_dof_with_counter_face(hashmap::unordered_map<tuple<T>, 
                                                          std::pair<emi_face<T,tuple<T>>, emi_face<T,tuple<T>>>> & line_face,
                                                          hashmap::unordered_map<triple<T>, 
                                                          std::pair<emi_face<T,triple<T>>, emi_face<T,triple<T>>>> & tri_face,
                                                          hashmap::unordered_map<quadruple<T>, 
                                                          std::pair<emi_face<T,quadruple<T>>, emi_face<T,quadruple<T>>>> & quad_face,
                                                          hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, 
                                                          mesh_int_t> & map_vertex_tag_to_dof)
{
  MPI_Comm comm = SF_COMM;
  int      size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  // look for counter part line_inteface to add to the map key<vertex,tag> -> value<dof>
  for(const auto & v : line_face) {
 
    // if the ranks are equal, then it has already added to the map
    if(v.second.first.rank == v.second.second.rank){
      continue;
    }
    // otherwise we look for the counter part with different rank
    emi_face<T,tuple<T>> surf_neighbor = 
      (v.second.first.rank != rank) ? v.second.first : v.second.second;
 
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(surf_neighbor.points.v1,surf_neighbor.tag);
      map_vertex_tag_to_dof.insert({Index_tag_old,-1});
    }
 
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(surf_neighbor.points.v2,surf_neighbor.tag);
      map_vertex_tag_to_dof.insert({Index_tag_old,-1});
    }
  }
 
  // look for counter part line_inteface to add to the map key<vertex,tag> -> value<dof>
  for(const auto & v : tri_face) {
 
    // if the ranks are equal, then it has already added to the map
    if(v.second.first.rank == v.second.second.rank){
      continue;
    }
    // otherwise we look for the counter part with different rank
    emi_face<T,triple<T>> surf_neighbor = 
      (v.second.first.rank != rank) ? v.second.first : v.second.second;
 
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(surf_neighbor.points.v1,surf_neighbor.tag);
      map_vertex_tag_to_dof.insert({Index_tag_old,-1});
    }
 
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(surf_neighbor.points.v2,surf_neighbor.tag);
      map_vertex_tag_to_dof.insert({Index_tag_old,-1});
    }
 
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(surf_neighbor.points.v3,surf_neighbor.tag);
      map_vertex_tag_to_dof.insert({Index_tag_old,-1});
    }
  }
 
  // look for counter part line_inteface to add to the map key<vertex,tag> -> value<dof>
  for(const auto & v : quad_face) {
 
    // if the ranks are equal, then it has already added to the map
    if(v.second.first.rank == v.second.second.rank){
      continue;
    }
    // otherwise we look for the counter part with different rank
    emi_face<T,quadruple<T>> qsurf_neighbor = 
      (v.second.first.rank != rank) ? v.second.first : v.second.second;
 
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(qsurf_neighbor.points.v1,qsurf_neighbor.tag);
      map_vertex_tag_to_dof.insert({Index_tag_old,-1});
    }
 
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(qsurf_neighbor.points.v2,qsurf_neighbor.tag);
      map_vertex_tag_to_dof.insert({Index_tag_old,-1});
    }
 
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(qsurf_neighbor.points.v3,qsurf_neighbor.tag);
      map_vertex_tag_to_dof.insert({Index_tag_old,-1});
    }
 
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(qsurf_neighbor.points.v4,qsurf_neighbor.tag);
      map_vertex_tag_to_dof.insert({Index_tag_old,-1});
    }
  }
 
  // assign the counter faces which we haven't assigned it yet the new dof
  assign_dof_on_counter_face(map_vertex_tag_to_dof,comm);
}
 
template<class T, class S> inline
void compute_map_vertex_to_dof(meshdata<T,S> & mesh,
                                  const SF_nbr numbering,
                                  hashmap::unordered_map<T,T> & vertex2ptsdata,
                                  hashmap::unordered_set<int> & extra_tags,
                                  hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>,
                                                         mesh_int_t> & map_vertex_tag_to_dof)
{
  MPI_Comm comm = mesh.comm;
  int      size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  hashmap::unordered_map<mesh_int_t, intersection_data> map_vertex_to_tags_data_ranks;
 
  const T* con = mesh.con.data();
  const T* nbr = mesh.get_numbering(numbering).data();
  const SF::vector<T> & rnod = mesh.get_numbering(numbering);
 
  hashmap::unordered_map<T,T> g2ptsData;
 
  for(size_t i=0; i<mesh.con.size(); i++)
  {
    g2ptsData[rnod[con[i]]] = vertex2ptsdata[rnod[con[i]]];
  }
 
  //-----------------------------------------------------------------
  // Step 1: Local Information Gathering
  //-----------------------------------------------------------------
  //    - Identify Simple Cases: For any vertex that is not on an interface (i.e., it's completely inside the extracellular space or completely
  //      inside a unique myocyte), the mapping is simple. The new DOF index is the same as the original vertex index. The function populates the
  //      output map, map_vertex_tag_to_dof, with these direct mappings, e.g., (vertex {1}, tag_3) -> 1.
  //    - Collect Intersection Data: This is the most important part of this step. For every vertex, it builds a list of all the regions (tags) and
  //      MPI ranks that share it, based on the elements the current process knows about. This information is stored in a temporary map called
  //      map_vertex_to_tags_data_ranks.
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++)
  {
    T tag = mesh.tag[eidx];
    auto pos_extra = extra_tags.find(tag);
    for (int n = mesh.dsp[eidx]; n < mesh.dsp[eidx+1];n++)
    {
      T gIndex = rnod[con[n]];
      T data_on_gIndex = g2ptsData[rnod[con[n]]];
 
      // This handles inner DoFs, either within the extracellular domain or vertices belonging to the intracellular domain.
      // and marked as a visited dofs and saved key<gIndex,tag)> -> val<gIndex>
      if((pos_extra != extra_tags.end()) || data_on_gIndex==0){
        std::pair <mesh_int_t,mesh_int_t> Index_tag_old = std::make_pair(gIndex,tag);
        if (map_vertex_tag_to_dof.find(Index_tag_old) == map_vertex_tag_to_dof.end() ) 
        {
          map_vertex_tag_to_dof.insert({Index_tag_old,gIndex});
        }
      }
 
      // When gIndex is located on an interface — membrane, gap junction, or their intersection.
      auto it = map_vertex_to_tags_data_ranks.find(gIndex);
      if (it != map_vertex_to_tags_data_ranks.end() ) 
      {
        intersection_data& value = it->second;
        bool check_tag_rank = true;
        for(int i=0; i<MAX_INTERSECTIONS; ++i) {
          if(value.tags[i] == tag && value.ranks[i] == rank) {
            check_tag_rank = false;
            break;
          }
        }
 
        // Insert the tag and corresponding rank if they haven’t been inserted yet.
        if(check_tag_rank){
          for(int i=0; i<MAX_INTERSECTIONS; ++i) {
            if(value.tags[i] == -1 && value.ranks[i] == -1) {
              value.tags[i] = tag;
              value.data[i] = data_on_gIndex;
              value.ranks[i] = rank;
              break;
            }
          }
        }
      }
      else
      {
        // insert gIndex to the map_vertex_to_tags_data_ranks
        intersection_data value;
        value.tags[0] = tag;
        value.data[0] = data_on_gIndex;
        value.ranks[0] = rank;
        map_vertex_to_tags_data_ranks.insert({gIndex,value});
      }
    }
  }
  // until here: map_vertex_to_tags_data_ranks knows about the interfaces for the elements it owns, 
  // but it doesn't know if a vertex it owns is also part of an interface on another process.
 
  //-----------------------------------------------------------------
  // Step 2: Globalize Information with MPI
  //-----------------------------------------------------------------
  //    - assign_counter_vertices_tuple: This helper function is called to manage the communication. It takes the partial
  //      map_vertex_to_tags_data_ranks from each process and uses MPI (specifically MPI_Exchange, which is likely a wrapper around routines like
  //      MPI_Alltoallv) to send and receive data.
  //    - Build Global Map: After the communication, the partial maps are merged. The result is that the map_vertex_to_tags_data_ranks on every 
  //      process now contains the complete sharing information for every vertex in the entire mesh. Each process now knows the full list of tags
  //      and ranks associated with any given vertex.
  assign_counter_vertices_tuple(map_vertex_to_tags_data_ranks, comm);
 
  //-----------------------------------------------------------------
  // Step 3: Identify and Count DOFs to Be Created
  //-----------------------------------------------------------------
  // Now that every process has the same global information, they can independently and deterministically decide which new DOFs to create.
  //   - Analyze Intersections: The code iterates through the now-global map_vertex_to_tags_data_ranks. For each vertex, it inspects the list of
  //     tags that share it to classify the type of interface.
  //   - Apply Rules: Based on the interface type, it decides which side needs a new, unique DOF index.
  //     * Membrane (Myocyte-Extracellular): The extracellular DOF keeps the original vertex index. The myocyte DOF is marked as needing a new, unique index.
  //     * Gap Junction (Myocyte-Myocyte): To separate the two myocytes, one of them needs a new index. A deterministic rule (e.g., the myocyte
  //       with the higher tag number gets the new index) is applied to ensure all processes make the same decision.
  //     * Complex Junctions: For even more complex intersections (e.g., where multiple myocytes and the extracellular space meet), similar
  //       rules are applied to create the necessary new DOFs.
  //   - Count Local Contribution: Each process counts how many new DOFs (shift) it is responsible for creating based on these rules.
  T shift = 0;
  hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, bool>  map_mark_new_dofs;
  for(const auto & key : map_vertex_to_tags_data_ranks) 
  {
    T gIndex  = key.first;
    const intersection_data& value = key.second;
    T data_on_gIndex = g2ptsData[gIndex];
 
    // find the first -1 in value
    int first_index = 0;
    for (int i = 0; i < MAX_INTERSECTIONS; ++i) {
      if(value.tags[i] == -1 && value.ranks[i] == -1) {
        first_index = i;
        break;
      }else{
        // check if all the intersection_data has the same g2ptsData
        if(value.data[i]!=data_on_gIndex)
        {
          // Error for unhandled cases.
          fprintf(stderr,
          "WARN: Due to an inconsistency in ptsData computation at gIndex %lld with tag %lld, "
          "the program is stopped while preparing the data structure to introduce DOFs.\n",
          static_cast<long long>(gIndex), static_cast<long long>(value.tags[i]));
          fprintf(stderr, "\n");
          EXIT(1);
        }
      }
      if (i == MAX_INTERSECTIONS -1) first_index = MAX_INTERSECTIONS;
    }
 
    if(data_on_gIndex==1){ // membrane
      // check the mumber of tag numbers, 
      int count_intra_tags = 0;
      int count_extra_tags = 0;
      for (int i = 0; i < first_index; ++i) {
        if(value.tags[i]!=-1 && extra_tags.find(value.tags[i])==extra_tags.end()){
          count_intra_tags += 1;
        }else if(value.tags[i]!=-1 && extra_tags.find(value.tags[i])!=extra_tags.end()){
          count_extra_tags += 1;
        }
      }
 
      // Check if the collected data at gIndex is indeed on the membrane.
      // A valid membrane vertex must intersect exactly one intracellular tag and one or more extracellular tags.
      if(count_intra_tags>1 || count_extra_tags == 0)
      {
        // Error for unhandled cases.
        fprintf(stderr,
        "WARN: Due to an issue in defining ptsData for gIndex %lld, , on the membrane"
        "the program is stopped while preparing the data structure to introduce DOFs.\n",
        static_cast<long long>(gIndex));
        fprintf(stderr, "\n");
        EXIT(1);
      }
 
      int index_intra = -1;
      // Find the intracellular tag number. Since data_on_gIndex == 1, there will be only one tag corresponding to the intracellular domain.
      for (int i = 0; i < first_index; ++i) {
        if(value.tags[i]!=-1 && extra_tags.find(value.tags[i])==extra_tags.end()){
          index_intra = i;
          break;
        }
      }
      if(index_intra == -1) { 
        fprintf(stderr, "WARN: unhandled case in membrane while decoupling: gIndex %lld",
                static_cast<long long>(gIndex));
        fprintf(stderr, "\n");
        EXIT(1);
      }
      if(index_intra != -1) {
        T  tag_myocyte = value.tags[index_intra];
        T  rank_myocyte = value.ranks[index_intra];
 
        if(rank_myocyte == rank) { // mark intracellular tag with gIndex
          std::pair <mesh_int_t,mesh_int_t> Index_tag_next = std::make_pair(gIndex,tag_myocyte);
          if (map_mark_new_dofs.find(Index_tag_next) == map_mark_new_dofs.end()) {
            map_mark_new_dofs.insert({Index_tag_next,true}); //mark a new index
            shift++;
          }
        }
      }
    }
    else if(data_on_gIndex==2){ // gap junction
      // count the number of extra tags and intra tags
      int count_intra_tags = 0;
      int count_extra_tags = 0;
      for (int i = 0; i < first_index; ++i) {
        if(value.tags[i]!=-1 && extra_tags.find(value.tags[i])==extra_tags.end()){
          count_intra_tags += 1;
        }else if(value.tags[i]!=-1 && extra_tags.find(value.tags[i])!=extra_tags.end()){
          count_extra_tags += 1;
        }
      }
      if(count_intra_tags!=2 || count_extra_tags != 0){
        // Error for unhandled cases.
        fprintf(stderr,
        "WARN: Due to an issue in defining ptsData for gIndex %lld, on the gapjunction"
        "the program is stopped while preparing the data structure to introduce DOFs.\n",
        static_cast<long long>(gIndex));
        fprintf(stderr, "\n");
        EXIT(1);
      }
 
      T tag1 = -1;
      T tag2 = -1;
      T rank1 = -1;
      T rank2 = -1;
      bool tag1_checked = false;
      bool tag2_checked = false;
      // find the tag1 and tag2 on the gap junction
      for (int i = 0; i < first_index; ++i)
      {
        auto pos_extra = extra_tags.find(value.tags[i]);
        if(pos_extra == extra_tags.end()) {
          if(tag1==-1){
            tag1 = value.tags[i];
            rank1 = value.ranks[i];
            tag1_checked = true;
          }else if(tag1!=-1){
            tag2 = value.tags[i];
            rank2 = value.ranks[i];
            tag2_checked = true;
            break;
          }
        }
      }
 
      if(!tag1_checked && !tag2_checked || (rank1==rank2 && rank1!=rank) ) { 
        fprintf(stderr,
                "WARN: unhandled case in gap junction while decoupling: gIndex %lld with tags %lld:%lld ",
                static_cast<long long>(gIndex),
                static_cast<long long>(tag1),
                static_cast<long long>(tag2));
        fprintf(stderr, "\n");
        EXIT(1);
      }
 
      if(rank1==rank2 && (rank1 == rank)){ // If both tag numbers belong to the same rank, keep the smaller tag’s vertex index and reserve a new index for the larger tag.
          T smallest_tag = tag1 < tag2? tag1:tag2;
          std::pair <mesh_int_t,mesh_int_t> Index_tag_old = std::make_pair(gIndex,smallest_tag);
          if (map_vertex_tag_to_dof.find(Index_tag_old) == map_vertex_tag_to_dof.end() ) {
              map_vertex_tag_to_dof.insert({Index_tag_old,gIndex});
          }
 
          T bigger_tag = tag1 < tag2? tag2:tag1;
          std::pair <mesh_int_t,mesh_int_t> Index_tag_new = std::make_pair(gIndex,bigger_tag);
          if (map_mark_new_dofs.find(Index_tag_new) == map_mark_new_dofs.end() ) {
              map_mark_new_dofs.insert({Index_tag_new,true});
              shift++;
          }
      }
      else if(rank1!=rank2 && (rank1 == rank)) { // If one belong to the current rank, keep the smaller tag’s vertex index else reserve a new index for the larger tag.
          if(tag1<tag2) {
              std::pair <mesh_int_t,mesh_int_t> Index_tag_old = std::make_pair(gIndex,tag1);
              if (map_vertex_tag_to_dof.find(Index_tag_old) == map_vertex_tag_to_dof.end() ) {
                  map_vertex_tag_to_dof.insert({Index_tag_old,gIndex});
              }
          } else {
              std::pair <mesh_int_t,mesh_int_t> Index_tag_new = std::make_pair(gIndex,tag1);
              if (map_mark_new_dofs.find(Index_tag_new) == map_mark_new_dofs.end() ) {
                  map_mark_new_dofs.insert({Index_tag_new,true});
                  shift++;
              }
          }
      }
      else if(rank1!=rank2 && (rank2 == rank)) { // If one belong to the current rank, keep the smaller tag’s vertex index else reserve a new index for the larger tag.
          if(tag2<tag1) {
              std::pair <mesh_int_t,mesh_int_t> Index_tag_old = std::make_pair(gIndex,tag2);
              if (map_vertex_tag_to_dof.find(Index_tag_old) == map_vertex_tag_to_dof.end() ) {
                  map_vertex_tag_to_dof.insert({Index_tag_old,gIndex});
              }
          } else {
              std::pair <mesh_int_t,mesh_int_t> Index_tag_new = std::make_pair(gIndex,tag2);
              if (map_mark_new_dofs.find(Index_tag_new) == map_mark_new_dofs.end() ) {
                  map_mark_new_dofs.insert({Index_tag_new,true});
                  shift++;
              }
          }
      }
    }
    else if(data_on_gIndex==3)
    { 
      // count the number of extra tags and intra tags
      int count_intra_tags = 0;
      int count_extra_tags = 0;
      for (int i = 0; i < first_index; ++i) {
        if(value.tags[i]!=-1 && extra_tags.find(value.tags[i])==extra_tags.end()){
          count_intra_tags += 1;
        }else if(value.tags[i]!=-1 && extra_tags.find(value.tags[i])!=extra_tags.end()){
          count_extra_tags += 1;
        }
      }
      if(count_intra_tags!=2 || count_extra_tags == 0){
        // Error for unhandled cases.
        fprintf(stderr,
        "WARN: Due to an issue in defining ptsData for gIndex %lld, on intersection between membrane and gapjunction"
        "the program is stopped while preparing the data structure to introduce DOFs.\n",
        static_cast<long long>(gIndex));
        fprintf(stderr, "\n");
        EXIT(1);
      }
 
      // All extracellular tags have already been added in the previous loop that iterated over the mesh.
      // on the interface between gap juntion and membarne
      for (int i = 0; i < first_index; ++i) {
        // find the first two tag number which belong to itracellular tags, add for indexing
        auto pos_extra = extra_tags.find(value.tags[i]);
        if(value.tags[i]!=-1 && value.ranks[i] == rank && pos_extra == extra_tags.end()) { // myocytes, for sure we should have two myocytes
          std::pair <mesh_int_t,mesh_int_t> Index_tag_new = std::make_pair(gIndex,value.tags[i]);
          if (map_mark_new_dofs.find(Index_tag_new) == map_mark_new_dofs.end() )
          {
            map_mark_new_dofs.insert({Index_tag_new,true}); //mark a new index
            shift++;
          }
        }
      }
    }
  }
 
  vector<T> dsp_dof(size);
  MPI_Allgather(&shift, sizeof(T), MPI_BYTE, dsp_dof.data(), sizeof(T), MPI_BYTE, comm);
 
  //-----------------------------------------------------------------
  // Step 4: Assign Unique Global Indices in Parallel
  //-----------------------------------------------------------------
  // The final step is to assign the new indices without any conflicts between processes.
  //   - Calculate Global Offset (`start`): An MPI_Allgather is used to share the shift counts among all processes. Each process then calculates
  //     its own unique start index for the block of new DOFs it will create. This is done by taking the total number of vertices in the original
  //     mesh and adding the shift counts from all processes with a smaller rank. This is a parallel prefix sum, and it guarantees that the ranges
  //     of new indices created by different processes will not overlap.
  //   - Assign New Indices: Finally, the function loops through the vertices it marked for new DOF creation and assigns them a unique global
  //     index (newIndex = start + count). This final mapping, (original_vertex_index, tag) -> new_unique_dof_index, is inserted into the output
  //     map map_vertex_tag_to_dof.
  T start;
  if(rank==0){
    start = mesh.g_numpts;
  }
  else
  {
    start = mesh.g_numpts;
    for (int r = 0; r < rank; ++r)
    {
        start+= dsp_dof[r];
    }
  }
 
  T count = 0;
  auto lexicographic_comp_pair = [](const std::pair<mesh_int_t,mesh_int_t> & a, const std::pair<mesh_int_t,mesh_int_t> & b)
  {
    if (a.second == b.second) return a.first < b.first;
    return a.second < b.second;
  };
 
  map_mark_new_dofs.sort(lexicographic_comp_pair);
 
  for(const auto & entry : map_mark_new_dofs)
  {
    T newIndex = start+count;
    map_vertex_tag_to_dof.insert({entry.first,newIndex});
    count++;
  }
}
 
template<class K, class IntersectionIndices> inline
void assign_counter_dofs(hashmap::unordered_map<K, IntersectionIndices>& map, const MPI_Comm comm)
{
  size_t dsize = map.size();
  vector<K> key_vec(dsize);
  vector<IntersectionIndices> value_vec(dsize);
  IntersectionIndices indices;
  size_t idx = 0;
  for (const auto& v : map) {
    key_vec[idx] = v.first;
    value_vec[idx] = v.second;
    idx++;
  }
 
  vector<int> perm, dest;
  emi_select_merge_destinations(key_vec, dsize, sizeof(K) + sizeof(intersection_indices), dest, comm,
                                "assign_counter_dofs",
                                [](const K& key) { return hashmap::hash_ops<K>::hash(key); });
 
  commgraph<size_t> grph;
  grph.configure(dest, comm);
  size_t nrecv = sum(grph.rcnt);
 
  interval(perm, 0, dsize);
  binary_sort_copy(dest, perm);
 
  // fill send buffer and communicate
  vector<K> sbuff_key(dsize), rbuff_key(nrecv);
  vector<IntersectionIndices> sbuff_value(dsize), ibuff_value(dsize), rbuff_value(nrecv);
  for (size_t i = 0; i < dsize; i++) {
    sbuff_key[i] = key_vec[perm[i]];
    sbuff_value[i] = value_vec[perm[i]];
  }
  MPI_Exchange(grph, sbuff_key, rbuff_key, comm);
  MPI_Exchange(grph, sbuff_value, rbuff_value, comm);
 
  hashmap::unordered_map<K, IntersectionIndices> rmap;
  for (size_t i = 0; i < nrecv; i++) 
  {
    auto it = rmap.find(rbuff_key[i]);
    if (it != rmap.end()) 
    {
      IntersectionIndices& map_indices = it->second;
      IntersectionIndices& received_indices = rbuff_value[i];
 
      for (mesh_int_t received_index : received_indices.indices)
      {
        if (received_index == -1) continue;
        bool found = false;
        for (mesh_int_t map_index : map_indices.indices)
        {
          if (map_index == received_index) {
            found = true;
            break;
          }
        }
        if (!found) 
        {
          for (size_t k = 0; k < MAX_INTERSECTIONS; ++k) {
            if (map_indices.indices[k] == -1) {
              map_indices.indices[k] = received_index;
              break;
            }
          }
        }
      }
    }
    else
    {
      rmap.insert({ rbuff_key[i], rbuff_value[i] });
    }
  }
 
  for (size_t i = 0; i < nrecv; i++) {
    auto it = rmap.find(rbuff_key[i]);
    if (it != rmap.end()) rbuff_value[i] = it->second;
  }
 
  // sending the assigned data back to original rank
  grph.transpose();
  MPI_Exchange(grph, rbuff_value, ibuff_value, comm);
 
  for (size_t i = 0; i < dsize; i++) {
      auto it = map.find(sbuff_key[i]);
      if (it != map.end()) it->second = ibuff_value[i];
  }
}
template<class T, class S> inline
void update_emi_mesh_with_dofs(meshdata<T,S> & mesh, 
                                 const SF_nbr numbering,
                                 hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, 
                                                        mesh_int_t> & map_vertex_tag_to_dof,
                                 hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex)
{
  // map between old index to set of new indices
  hashmap::unordered_map<mesh_int_t,intersection_indices> vertex2dof;
  MPI_Comm comm = mesh.comm;
  mesh.globalize(numbering);
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    T tag = mesh.tag[eidx];
    T size_elem = mesh.dsp[eidx+1]-mesh.dsp[eidx];
    
    for (int n = mesh.dsp[eidx]; n < mesh.dsp[eidx+1];n++)
    {
      T gIndex = mesh.con[n];
 
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(gIndex,tag);
 
      auto it_new = map_vertex_tag_to_dof.find(Index_tag_old);
      if (it_new != map_vertex_tag_to_dof.end() ) 
      {
        mesh_int_t dof_new = (*it_new).second;
        dof2vertex.insert({dof_new, gIndex});
        mesh.con[n] = dof_new;
 
        auto it = vertex2dof.find(gIndex);
        if (it != vertex2dof.end())
        {
          intersection_indices& indices = it->second;
          bool found = false;
          for(T t : indices.indices) {
            if(t == dof_new) {
              found = true;
              break;
            }
          }
          if(!found) {
            // the default size is 3, since the maximum intersection shoudl be on the intersection between the gap junction and memeranre, and since the 2 indiecs for myocyesand one indices for extracellular.
            for(size_t i=0; i<MAX_INTERSECTIONS; ++i) {  
              if(indices.indices[i] == -1) {
                indices.indices[i] = dof_new;
                break;
              }
            }
          }
        }else{
          intersection_indices indices;
          indices.indices[0] = dof_new;
          vertex2dof.insert({gIndex,indices});
        }
      }
    }
  }
  
  assign_counter_dofs(vertex2dof, comm);
 
  // Now update dof2vertex using vertex2dof
  for (const auto& [old_idx, indices] : vertex2dof) {
    for (mesh_int_t new_idx : indices.indices) {
      if (new_idx != -1) {
        // Avoid overwriting unless necessary
        if (dof2vertex.find(new_idx) == dof2vertex.end()) {
          dof2vertex.insert({new_idx, old_idx});
        }
      }
    }
  }
 
  mesh.localize(numbering);
}
 
template<class T, class S> inline
void update_map_indices_to_petsc(meshdata<T,S> & mesh, 
                                 const SF_nbr numbering_ref,const SF_nbr numbering_petsc,
                                 const hashmap::unordered_set<int> & extra_tags,
                                 hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, 
                                                         std::pair<mesh_int_t,mesh_int_t>> & map_vertex_tag_to_dof_petsc,
                                 hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex,
                                 SF::vector<mesh_int_t> & elemTag_emi_mesh)
{
  const SF::vector<T> & ref_nbr = mesh.get_numbering(numbering_ref);
  const SF::vector<T> & petsc_nbr = mesh.get_numbering(numbering_petsc);
  
  elemTag_emi_mesh.resize(mesh.l_numelem);
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    T tag = mesh.tag[eidx];
    elemTag_emi_mesh[eidx] = 1; // assigned 1 for for extracellular located on extracellular side
    if(extra_tags.find(tag) == extra_tags.end())
      elemTag_emi_mesh[eidx] = 2; // assigned 2 for for intracellular located on myocyte 
    for (int n = mesh.dsp[eidx]; n < mesh.dsp[eidx+1];n++)
    {
      T l_Indx = mesh.con[n];
      T petsc_Idx = petsc_nbr[l_Indx];
      T n_Indx = ref_nbr[l_Indx];
      T o_Indx = dof2vertex[n_Indx];
 
      std::pair <mesh_int_t,mesh_int_t> oIndx_tag;
      oIndx_tag = std::make_pair(o_Indx,tag);
      
      auto it_new = map_vertex_tag_to_dof_petsc.find(oIndx_tag);
      if (it_new != map_vertex_tag_to_dof_petsc.end() ) 
      {
        std::pair <mesh_int_t,mesh_int_t> nIndx_petsc;
        nIndx_petsc = std::make_pair(n_Indx,petsc_Idx);
 
        (*it_new).second = nIndx_petsc;
      }
    }
  }
}
 
template<class T, class S>
inline void update_faces_on_surface_mesh_after_decoupling_with_dofs(meshdata<T, S> & mesh, 
                                            hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, 
                                                                   mesh_int_t> & map_vertex_tag_to_dof,
                                            hashmap::unordered_map<tuple<T>, 
                                                                   std::pair<emi_face<T,tuple<T>>, emi_face<T,tuple<T>>>> & line_face,
                                            hashmap::unordered_map<triple<T>, 
                                                                   std::pair<emi_face<T,triple<T>>, emi_face<T,triple<T>>>> & tri_face,
                                            hashmap::unordered_map<quadruple<T>, 
                                                                   std::pair<emi_face<T,quadruple<T>>, emi_face<T,quadruple<T>>>> & quad_face)
{
  MPI_Comm comm = SF_COMM;
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  SF::vector<mesh_int_t> idxbuff(mesh.con);
  binary_sort(idxbuff); unique_resize(idxbuff);
 
  hashmap::unordered_map<T,T> g2l;
  hashmap::unordered_map<T,T> l2g;
  for(size_t i=0; i<idxbuff.size(); i++){
    g2l[idxbuff[i]] = i;
    l2g[i] = idxbuff[i];
  }
 
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    T tag = mesh.tag[eidx];
    T size_elem = mesh.dsp[eidx+1]-mesh.dsp[eidx];
    std::vector<mesh_int_t> elem_nodes;
    for (int n = mesh.dsp[eidx]; n < mesh.dsp[eidx+1];n++)
    {
      T gIndex = mesh.con[n];
      elem_nodes.push_back(gIndex);
    }
    std::sort(elem_nodes.begin(),elem_nodes.end());
    if(elem_nodes.size()==2){
      tuple<T> key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      std::pair<emi_face<T,tuple<T>>, emi_face<T,tuple<T>>> value = line_face[key];
 
      line_face.erase(key);
      std::vector<mesh_int_t> new_nodes(2);
 
      auto tag_key = (value.first.rank == rank) ? value.first.tag : value.second.tag;
      {
        std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
        Index_tag_old = std::make_pair(elem_nodes[0],tag_key);
        mesh_int_t Index_new = map_vertex_tag_to_dof[Index_tag_old];
        new_nodes[0] = Index_new;
      }
      {
        std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
        Index_tag_old = std::make_pair(elem_nodes[1],tag_key);
        mesh_int_t Index_new = map_vertex_tag_to_dof[Index_tag_old];
        new_nodes[1] = Index_new;
      }
      std::sort(new_nodes.begin(),new_nodes.end());
      // update key
      tuple<T> new_key;
      new_key.v1 = new_nodes[0];
      new_key.v2 = new_nodes[1];
      
      // first
      {
        mesh_int_t Index_new1 = map_vertex_tag_to_dof[std::make_pair(value.first.points.v1, value.first.tag)];
        mesh_int_t Index_new2 = map_vertex_tag_to_dof[std::make_pair(value.first.points.v2, value.first.tag)];
        value.first.points.v1 = Index_new1; 
        value.first.points.v2 = Index_new2; 
      }
      // second
      {
        mesh_int_t Index_new1 = map_vertex_tag_to_dof[std::make_pair(value.second.points.v1, value.second.tag)];
        mesh_int_t Index_new2 = map_vertex_tag_to_dof[std::make_pair(value.second.points.v2, value.second.tag)];
        value.second.points.v1 = Index_new1; 
        value.second.points.v2 = Index_new2; 
      }
      line_face.insert({new_key,value});
    }
    if(elem_nodes.size()==3){
      triple<T> key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      std::pair<emi_face<T,triple<T>>, emi_face<T,triple<T>>> value = tri_face[key];
      tri_face.erase(key);
      std::vector<mesh_int_t> new_nodes(3);
 
      auto tag_key = (value.first.rank == rank) ? value.first.tag : value.second.tag;
      {
        std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
        Index_tag_old = std::make_pair(elem_nodes[0],tag_key);
        mesh_int_t Index_new = map_vertex_tag_to_dof[Index_tag_old];
        new_nodes[0] = Index_new;
      }
      {
        std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
        Index_tag_old = std::make_pair(elem_nodes[1],tag_key);
        mesh_int_t Index_new = map_vertex_tag_to_dof[Index_tag_old];
        new_nodes[1] = Index_new;
      }
      {
        std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
        Index_tag_old = std::make_pair(elem_nodes[2],tag_key);
        mesh_int_t Index_new = map_vertex_tag_to_dof[Index_tag_old];
        new_nodes[2] = Index_new;
      }
      std::sort(new_nodes.begin(),new_nodes.end());
      // update key
      triple<T> new_key;
      new_key.v1 = new_nodes[0];
      new_key.v2 = new_nodes[1];
      new_key.v3 = new_nodes[2];
      
      // first
      {
        mesh_int_t Index_new1 = map_vertex_tag_to_dof[std::make_pair(value.first.points.v1, value.first.tag)];
        mesh_int_t Index_new2 = map_vertex_tag_to_dof[std::make_pair(value.first.points.v2, value.first.tag)];
        mesh_int_t Index_new3 = map_vertex_tag_to_dof[std::make_pair(value.first.points.v3, value.first.tag)];
        value.first.points.v1 = Index_new1; 
        value.first.points.v2 = Index_new2;
        value.first.points.v3 = Index_new3; 
      }
      // second
      {
        mesh_int_t Index_new1 = map_vertex_tag_to_dof[std::make_pair(value.second.points.v1, value.second.tag)];
        mesh_int_t Index_new2 = map_vertex_tag_to_dof[std::make_pair(value.second.points.v2, value.second.tag)];
        mesh_int_t Index_new3 = map_vertex_tag_to_dof[std::make_pair(value.second.points.v3, value.second.tag)];
        value.second.points.v1 = Index_new1; 
        value.second.points.v2 = Index_new2;
        value.second.points.v3 = Index_new3; 
      }
      tri_face.insert({new_key,value});
    }
    if(elem_nodes.size()==4){
      quadruple<T> key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      key.v4 = elem_nodes[3];
      std::pair<emi_face<T,quadruple<T>>, emi_face<T,quadruple<T>>> value = quad_face[key];
      quad_face.erase(key);
      std::vector<mesh_int_t> new_nodes(4);
 
      auto tag_key = (value.first.rank == rank) ? value.first.tag : value.second.tag;
      {
        std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
        Index_tag_old = std::make_pair(elem_nodes[0],tag_key);
        mesh_int_t Index_new = map_vertex_tag_to_dof[Index_tag_old];
        new_nodes[0] = Index_new;
      }
      {
        std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
        Index_tag_old = std::make_pair(elem_nodes[1],tag_key);
        mesh_int_t Index_new = map_vertex_tag_to_dof[Index_tag_old];
        new_nodes[1] = Index_new;
      }
      {
        std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
        Index_tag_old = std::make_pair(elem_nodes[2],tag_key);
        mesh_int_t Index_new = map_vertex_tag_to_dof[Index_tag_old];
        new_nodes[2] = Index_new;
      }
      {
        std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
        Index_tag_old = std::make_pair(elem_nodes[3],tag_key);
        mesh_int_t Index_new = map_vertex_tag_to_dof[Index_tag_old];
        new_nodes[3] = Index_new;
      }
      std::sort(new_nodes.begin(),new_nodes.end());
      // update key
      quadruple<T> new_key;
      new_key.v1 = new_nodes[0];
      new_key.v2 = new_nodes[1];
      new_key.v3 = new_nodes[2];
      new_key.v4 = new_nodes[3];
 
      // first
      {
        mesh_int_t Index_new1 = map_vertex_tag_to_dof[std::make_pair(value.first.points.v1, value.first.tag)];
        mesh_int_t Index_new2 = map_vertex_tag_to_dof[std::make_pair(value.first.points.v2, value.first.tag)];
        mesh_int_t Index_new3 = map_vertex_tag_to_dof[std::make_pair(value.first.points.v3, value.first.tag)];
        mesh_int_t Index_new4 = map_vertex_tag_to_dof[std::make_pair(value.first.points.v4, value.first.tag)];
        value.first.points.v1 = Index_new1; 
        value.first.points.v2 = Index_new2;
        value.first.points.v3 = Index_new3; 
        value.first.points.v4 = Index_new4; 
      }
      // second
      {
        mesh_int_t Index_new1 = map_vertex_tag_to_dof[std::make_pair(value.second.points.v1, value.second.tag)];
        mesh_int_t Index_new2 = map_vertex_tag_to_dof[std::make_pair(value.second.points.v2, value.second.tag)];
        mesh_int_t Index_new3 = map_vertex_tag_to_dof[std::make_pair(value.second.points.v3, value.second.tag)];
        mesh_int_t Index_new4 = map_vertex_tag_to_dof[std::make_pair(value.second.points.v4, value.second.tag)];
        value.second.points.v1 = Index_new1; 
        value.second.points.v2 = Index_new2;
        value.second.points.v3 = Index_new3; 
        value.second.points.v4 = Index_new4; 
      }
      quad_face.insert({new_key,value});
    }
  }
}
 
template<class T, class S> inline
void compute_surface_mesh_with_counter_face(meshdata<T, S> & mesh, const SF_nbr numbering,
                                  hashmap::unordered_map<tuple<T>, 
                                  std::pair<emi_face<T,tuple<T>>, emi_face<T,tuple<T>>>> & line_face,
                                  hashmap::unordered_map<triple<T>, 
                                  std::pair<emi_face<T,triple<T>>, emi_face<T,triple<T>>>> & tri_face,
                                  hashmap::unordered_map<quadruple<T>, 
                                  std::pair<emi_face<T,quadruple<T>>, emi_face<T,quadruple<T>>>> & quad_face)
{
  mesh.register_numbering(numbering); 
 
  MPI_Comm comm = SF_COMM;
  int      size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  // To ensure deterministic order, we must sort the keys of the hashmaps before iterating.
  vector<tuple<T>> line_keys;
  for(auto const& [key, val] : line_face) line_keys.push_back(key);
  std::sort(line_keys.begin(), line_keys.end());
 
  vector<triple<T>> tri_keys;
  for(auto const& [key, val] : tri_face) tri_keys.push_back(key);
  std::sort(tri_keys.begin(), tri_keys.end());
 
  vector<quadruple<T>> quad_keys;
  for(auto const& [key, val] : quad_face) quad_keys.push_back(key);
  std::sort(quad_keys.begin(), quad_keys.end());
 
  vector<T> cnt(mesh.l_numelem);
  size_t idx = 0, cidx = 0;
 
  // Add faces to the surface mesh used in ionic computation.
  // For each {line,tri,quad}_face, add faces only if one of the pair matches the current rank,
  // or if both faces have the same rank, add both.
  for(const auto & key : line_keys) {
    const auto & v = line_face.at(key);
    bool both_faces = (v.first.rank == v.second.rank) ? true: false;
 
    cnt[idx] = 2;
 
    emi_face<T,tuple<T>> surf_first;
    emi_face<T,tuple<T>> surf_second;
    if(v.first.rank == rank)
    {
      surf_first = v.first;
      surf_second = v.second;
    } 
    else
    {
      surf_first = v.second;
      surf_second = v.first;       
    } 
 
    mesh.type[idx]     = Line;
    mesh.tag[idx]      = surf_first.tag;
    mesh.con[cidx + 0] = surf_first.points.v1;
    mesh.con[cidx + 1] = surf_first.points.v2;
    idx  += 1;
    cidx += 2;
 
    if(both_faces){
 
      cnt[idx] = 2;
      mesh.type[idx]     = Line;
      mesh.tag[idx]      = surf_second.tag;
      mesh.con[cidx + 0] = surf_second.points.v1;
      mesh.con[cidx + 1] = surf_second.points.v2;
 
      idx  += 1;
      cidx += 2;
    }
  }
 
  for(const auto & key : tri_keys) {
    const auto & v = tri_face.at(key);
    bool both_faces = (v.first.rank == v.second.rank) ? true: false;
 
    cnt[idx] = 3;
 
    emi_face<T,triple<T>> surf_first;
    emi_face<T,triple<T>> surf_second;
    if(v.first.rank == rank)
    {
      surf_first = v.first;
      surf_second = v.second;
    } 
    else
    {
      surf_first = v.second;
      surf_second = v.first;       
    } 
 
    mesh.type[idx]     = Tri;
    mesh.tag[idx]      = surf_first.tag;
    mesh.con[cidx + 0] = surf_first.points.v1;
    mesh.con[cidx + 1] = surf_first.points.v2;
    mesh.con[cidx + 2] = surf_first.points.v3;
 
    idx  += 1;
    cidx += 3;
 
    if(both_faces){
 
      cnt[idx] = 3;
 
      mesh.type[idx]     = Tri;
      mesh.tag[idx]      = surf_second.tag;
      mesh.con[cidx + 0] = surf_second.points.v1;
      mesh.con[cidx + 1] = surf_second.points.v2;
      mesh.con[cidx + 2] = surf_second.points.v3;
 
      idx  += 1;
      cidx += 3;
    }
  }
 
  for(const auto & key : quad_keys) {
    const auto & v = quad_face.at(key);
    bool both_faces = (v.first.rank == v.second.rank) ? true: false;
    cnt[idx] = 4;
 
    emi_face<T,quadruple<T>> surf_first;
    emi_face<T,quadruple<T>> surf_second;
    if(v.first.rank == rank)
    {
      surf_first = v.first;
      surf_second = v.second;
    } 
    else
    {
      surf_first = v.second;
      surf_second = v.first;     
    } 
 
    mesh.type[idx]     = Quad;
    mesh.tag[idx]      = surf_first.tag;
    mesh.con[cidx + 0] = surf_first.points.v1;
    mesh.con[cidx + 1] = surf_first.points.v2;
    mesh.con[cidx + 2] = surf_first.points.v3;
    mesh.con[cidx + 3] = surf_first.points.v4;
 
    idx  += 1;
    cidx += 4;
 
    if(both_faces){
 
      cnt[idx] = 4;
 
      mesh.type[idx]     = Quad;
      mesh.tag[idx]      = surf_second.tag;
      mesh.con[cidx + 0] = surf_second.points.v1;
      mesh.con[cidx + 1] = surf_second.points.v2;
      mesh.con[cidx + 2] = surf_second.points.v3;
      mesh.con[cidx + 3] = surf_second.points.v4;
 
      idx  += 1;
      cidx += 4;
    }
  }
  dsp_from_cnt(cnt, mesh.dsp);
}
 
template<class T, class S> inline
void compute_surface_mesh_with_unique_face(meshdata<T, S> & mesh, const SF_nbr numbering,
                                  hashmap::unordered_map<tuple<T>, 
                                  std::pair<emi_face<T,tuple<T>>, emi_face<T,tuple<T>>>> & line_face,
                                  hashmap::unordered_map<triple<T>, 
                                  std::pair<emi_face<T,triple<T>>, emi_face<T,triple<T>>>> & tri_face,
                                  hashmap::unordered_map<quadruple<T>, 
                                  std::pair<emi_face<T,quadruple<T>>, emi_face<T,quadruple<T>>>> & quad_face,
                                  hashmap::unordered_map<mesh_int_t, std::pair<mesh_int_t, mesh_int_t>> & map_elem_uniqueFace_to_tags)
{
  mesh.register_numbering(numbering); 
 
  MPI_Comm comm = SF_COMM;
  int      size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  // To ensure deterministic order, we must sort the keys of the hashmaps before iterating.
  vector<tuple<T>> line_keys;
  for(auto const& [key, val] : line_face) line_keys.push_back(key);
  std::sort(line_keys.begin(), line_keys.end());
 
  vector<triple<T>> tri_keys;
  for(auto const& [key, val] : tri_face) tri_keys.push_back(key);
  std::sort(tri_keys.begin(), tri_keys.end());
 
  vector<quadruple<T>> quad_keys;
  for(auto const& [key, val] : quad_face) quad_keys.push_back(key);
  std::sort(quad_keys.begin(), quad_keys.end());
 
  vector<T> cnt(mesh.l_numelem);
  size_t idx = 0, cidx = 0;
 
  // Add faces to the surface mesh used in ionic computation.
  // For each {line,tri,quad}_face, add faces only if one of the pair matches the current rank,
  // or if both faces have the same rank, add both.
  for(const auto & key : line_keys) {
    const auto & v = line_face.at(key);
    
    bool to_take_face = false;
    
    emi_face<T,tuple<T>> surf_take;
    std::pair <mesh_int_t,mesh_int_t> tag_pairs_take;
    {
      if(v.first.rank == rank && v.first.mark_to_take == true)
      {
        surf_take = v.first;
        to_take_face = true;
        tag_pairs_take = std::make_pair(v.first.tag,v.second.tag);
      } 
      else if(v.second.rank == rank && v.second.mark_to_take == true)
      {
        surf_take = v.second;
        to_take_face = true;
        tag_pairs_take = std::make_pair(v.second.tag,v.first.tag);        
      } 
 
      if(to_take_face)
      { 
        map_elem_uniqueFace_to_tags.insert({idx,tag_pairs_take});
        cnt[idx] = 2;
        mesh.type[idx]     = Line;
        mesh.tag[idx]      = surf_take.tag;
        mesh.con[cidx + 0] = surf_take.points.v1;
        mesh.con[cidx + 1] = surf_take.points.v2;
        idx  += 1;
        cidx += 2;
      } 
    }
  }
 
  for(const auto & key : tri_keys) {
    const auto & v = tri_face.at(key);
 
    bool to_take_face = false;
 
    emi_face<T,triple<T>> surf_take;
    std::pair <mesh_int_t,mesh_int_t> tag_pairs_take;
    {
      if(v.first.rank == rank && v.first.mark_to_take == true)
      {
        surf_take = v.first;
        to_take_face = true;
        tag_pairs_take = std::make_pair(v.first.tag,v.second.tag);
      } 
      else if(v.second.rank == rank && v.second.mark_to_take == true)
      {
        surf_take = v.second; 
        to_take_face = true;   
        tag_pairs_take = std::make_pair(v.second.tag,v.first.tag);  
      } 
 
      if(to_take_face)
      {
        map_elem_uniqueFace_to_tags.insert({idx,tag_pairs_take});
        cnt[idx] = 3;
        mesh.type[idx]     = Tri;
        mesh.tag[idx]      = surf_take.tag;
        mesh.con[cidx + 0] = surf_take.points.v1;
        mesh.con[cidx + 1] = surf_take.points.v2;
        mesh.con[cidx + 2] = surf_take.points.v3;
 
        idx  += 1;
        cidx += 3;
      }
    }
  }
 
  for(const auto & key : quad_keys) {
    const auto & v = quad_face.at(key);
 
    bool to_take_face = false;
 
    emi_face<T,quadruple<T>> surf_take;
    std::pair <mesh_int_t,mesh_int_t> tag_pairs_take;
    {
      if(v.first.rank == rank && v.first.mark_to_take == true)
      {
        surf_take = v.first;
        to_take_face = true;
        tag_pairs_take = std::make_pair(v.first.tag,v.second.tag);
      } 
      else if(v.second.rank == rank && v.second.mark_to_take == true)
      {
        surf_take = v.second; 
        to_take_face = true; 
        tag_pairs_take = std::make_pair(v.second.tag,v.first.tag);        
      } 
 
      if(to_take_face)
      {
        map_elem_uniqueFace_to_tags.insert({idx,tag_pairs_take});
        cnt[idx] = 4;
        mesh.type[idx]     = Quad;
        mesh.tag[idx]      = surf_take.tag;
        mesh.con[cidx + 0] = surf_take.points.v1;
        mesh.con[cidx + 1] = surf_take.points.v2;
        mesh.con[cidx + 2] = surf_take.points.v3;
        mesh.con[cidx + 3] = surf_take.points.v4;
 
        idx  += 1;
        cidx += 4;
      }
    }
  }
  dsp_from_cnt(cnt, mesh.dsp);
}
 
template<class T> inline
void create_reverse_elem_mapping_between_surface_meshes(
    hashmap::unordered_map<tuple<T>, std::pair<emi_face<T,tuple<T>>, emi_face<T,tuple<T>>>>& line_face,
    hashmap::unordered_map<triple<T>, std::pair<emi_face<T,triple<T>>, emi_face<T,triple<T>>>>& tri_face,
    hashmap::unordered_map<quadruple<T>, std::pair<emi_face<T,quadruple<T>>, emi_face<T,quadruple<T>>>>& quad_face,
    SF::vector<T>& vec_one_to_both_face,
    MPI_Comm comm)
{
  int rank;
  MPI_Comm_rank(comm, &rank);
 
  // To ensure deterministic order, we must sort the keys of the hashmaps before iterating.
  vector<tuple<T>> line_keys;
  for(auto const& [key, val] : line_face) line_keys.push_back(key);
  std::sort(line_keys.begin(), line_keys.end());
 
  vector<triple<T>> tri_keys;
  for(auto const& [key, val] : tri_face) tri_keys.push_back(key);
  std::sort(tri_keys.begin(), tri_keys.end());
 
  vector<quadruple<T>> quad_keys;
  for(auto const& [key, val] : quad_face) quad_keys.push_back(key);
  std::sort(quad_keys.begin(), quad_keys.end());
 
  size_t surf_elem_idx = 0;
  size_t w_counter_elem_idx = 0;
 
  size_t lsize_line = 0;
  size_t lsize_tri = 0;
  size_t lsize_quad = 0;
 
  // Count how many both-face elements are local on this rank.
  // If both faces are on the same rank, the both-face mesh has two entries.
  for(const auto& key : line_keys) {
    const auto& v = line_face.at(key);
    const bool local_involved = (v.first.rank == rank) || (v.second.rank == rank);
    if (!local_involved) {
      continue;
    }
    if(v.first.rank == v.second.rank) lsize_line+=2;
    else lsize_line+=1;
  }
  for(const auto& key : tri_keys) {
    const auto& v = tri_face.at(key);
    const bool local_involved = (v.first.rank == rank) || (v.second.rank == rank);
    if (!local_involved) {
      continue;
    }
    if(v.first.rank == v.second.rank) lsize_tri+=2;
    else lsize_tri+=1;
  }
  for(const auto& key : quad_keys) {
    const auto& v = quad_face.at(key);
    const bool local_involved = (v.first.rank == rank) || (v.second.rank == rank);
    if (!local_involved) {
      continue;
    }
    if(v.first.rank == v.second.rank) lsize_quad+=2;
    else lsize_quad+=1;
  }
 
  vec_one_to_both_face.resize(lsize_line + lsize_tri + lsize_quad);
 
  // Fill mapping in deterministic order:
  // each both-face element gets the index of its corresponding one-face element.
  for (const auto& key : line_keys) {
    const auto& v = line_face.at(key);
    const bool local_involved = (v.first.rank == rank) || (v.second.rank == rank);
    if (!local_involved) {
      continue;
    }
    vec_one_to_both_face[w_counter_elem_idx++] = surf_elem_idx;
    if (v.first.rank == v.second.rank) {
        vec_one_to_both_face[w_counter_elem_idx++] = surf_elem_idx;
    }
    surf_elem_idx++;
  }
 
  for (const auto& key : tri_keys) {
    const auto& v = tri_face.at(key);
    const bool local_involved = (v.first.rank == rank) || (v.second.rank == rank);
    if (!local_involved) {
      continue;
    }
    vec_one_to_both_face[w_counter_elem_idx++] = surf_elem_idx;
    if (v.first.rank == v.second.rank) {
        vec_one_to_both_face[w_counter_elem_idx++] = surf_elem_idx;
    }
    surf_elem_idx++;
  }
 
  for (const auto& key : quad_keys) {
    const auto& v = quad_face.at(key);
    const bool local_involved = (v.first.rank == rank) || (v.second.rank == rank);
    if (!local_involved) {
      continue;
    }
    vec_one_to_both_face[w_counter_elem_idx++] = surf_elem_idx;
    if (v.first.rank == v.second.rank) {
        vec_one_to_both_face[w_counter_elem_idx++] = surf_elem_idx;
    }
    surf_elem_idx++;
  }
}
 
template<class T> inline 
void added_counter_faces_to_map(hashmap::unordered_map<tuple<T>, 
                                  std::pair<emi_face<T,tuple<T>>, emi_face<T,tuple<T>>>> & line_face,
                                  hashmap::unordered_map<triple<T>, 
                                  std::pair<emi_face<T,triple<T>>, emi_face<T,triple<T>>>> & tri_face,
                                  hashmap::unordered_map<quadruple<T>, 
                                  std::pair<emi_face<T,quadruple<T>>, emi_face<T,quadruple<T>>>> & quad_face)
{
  // Buffer inserts not to mutate an unordered_map while iterating over it and later add to unordered_map 
  std::vector<std::pair<tuple<T>, std::pair<emi_face<T,tuple<T>>, emi_face<T,tuple<T>>>>> line_inserts;
  for(const auto & v : line_face) {
    SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>> surf_first = v.second.first;
    SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>> surf_second = v.second.second;
 
    std::vector<mesh_int_t> elem_nodes_first; elem_nodes_first.resize(2);
    std::vector<mesh_int_t> elem_nodes_second; elem_nodes_second.resize(2);
 
    SF::tuple<mesh_int_t> key_first;
    elem_nodes_first[0] = surf_first.points.v1;
    elem_nodes_first[1] = surf_first.points.v2;
    std::sort(elem_nodes_first.begin(),elem_nodes_first.end());
    key_first.v1 = elem_nodes_first[0];
    key_first.v2 = elem_nodes_first[1];
 
    SF::tuple<mesh_int_t> key_second;
    elem_nodes_second[0] = surf_second.points.v1;
    elem_nodes_second[1] = surf_second.points.v2;
    std::sort(elem_nodes_second.begin(),elem_nodes_second.end());
    key_second.v1 = elem_nodes_second[0];
    key_second.v2 = elem_nodes_second[1];
 
    const bool same_key_first = (v.first.v1 == key_first.v1 && v.first.v2 == key_first.v2);
    if (!same_key_first && line_face.find(key_first) == line_face.end()) {
      line_inserts.push_back({key_first, v.second});
    }
    const bool same_key_second = (v.first.v1 == key_second.v1 && v.first.v2 == key_second.v2);
    if (!same_key_second && line_face.find(key_second) == line_face.end()) {
      line_inserts.push_back({key_second, v.second});
    }
  }
  for (const auto & entry : line_inserts) line_face.insert(entry);
 
  std::vector<std::pair<triple<T>, std::pair<emi_face<T,triple<T>>, emi_face<T,triple<T>>>>> tri_inserts;
  for(const auto & v : tri_face) {
    SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>> surf_first = v.second.first;
    SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>> surf_second = v.second.second;
 
    std::vector<mesh_int_t> elem_nodes_first; elem_nodes_first.resize(3);
    std::vector<mesh_int_t> elem_nodes_second; elem_nodes_second.resize(3);
 
    SF::triple<mesh_int_t> key_first;
    elem_nodes_first[0] = surf_first.points.v1;
    elem_nodes_first[1] = surf_first.points.v2;
    elem_nodes_first[2] = surf_first.points.v3;
    std::sort(elem_nodes_first.begin(),elem_nodes_first.end());
    key_first.v1 = elem_nodes_first[0];
    key_first.v2 = elem_nodes_first[1];
    key_first.v3 = elem_nodes_first[2];
 
    SF::triple<mesh_int_t> key_second;
    elem_nodes_second[0] = surf_second.points.v1;
    elem_nodes_second[1] = surf_second.points.v2;
    elem_nodes_second[2] = surf_second.points.v3;
    std::sort(elem_nodes_second.begin(),elem_nodes_second.end());
    key_second.v1 = elem_nodes_second[0];
    key_second.v2 = elem_nodes_second[1];
    key_second.v3 = elem_nodes_second[2];
 
    const bool same_key_first = (v.first.v1 == key_first.v1 &&
                                 v.first.v2 == key_first.v2 &&
                                 v.first.v3 == key_first.v3);
    if (!same_key_first && tri_face.find(key_first) == tri_face.end()) {
      tri_inserts.push_back({key_first, v.second});
    }
    const bool same_key_second = (v.first.v1 == key_second.v1 &&
                                  v.first.v2 == key_second.v2 &&
                                  v.first.v3 == key_second.v3);
    if (!same_key_second && tri_face.find(key_second) == tri_face.end()) {
      tri_inserts.push_back({key_second, v.second});
    }
  }
  for (const auto & entry : tri_inserts) tri_face.insert(entry);
 
  std::vector<std::pair<quadruple<T>, std::pair<emi_face<T,quadruple<T>>, emi_face<T,quadruple<T>>>>> quad_inserts;
  for(const auto & v : quad_face) {
    SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>> surf_first = v.second.first;
    SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>> surf_second = v.second.second;
 
    std::vector<mesh_int_t> elem_nodes_first; elem_nodes_first.resize(4);
    std::vector<mesh_int_t> elem_nodes_second; elem_nodes_second.resize(4);
 
    SF::quadruple<mesh_int_t> key_first;
    elem_nodes_first[0] = surf_first.points.v1;
    elem_nodes_first[1] = surf_first.points.v2;
    elem_nodes_first[2] = surf_first.points.v3;
    elem_nodes_first[3] = surf_first.points.v4;
    std::sort(elem_nodes_first.begin(),elem_nodes_first.end());
    key_first.v1 = elem_nodes_first[0];
    key_first.v2 = elem_nodes_first[1];
    key_first.v3 = elem_nodes_first[2];
    key_first.v4 = elem_nodes_first[3];
 
    SF::quadruple<mesh_int_t> key_second;
    elem_nodes_second[0] = surf_second.points.v1;
    elem_nodes_second[1] = surf_second.points.v2;
    elem_nodes_second[2] = surf_second.points.v3;
    elem_nodes_second[3] = surf_second.points.v4;
    std::sort(elem_nodes_second.begin(),elem_nodes_second.end());
    key_second.v1 = elem_nodes_second[0];
    key_second.v2 = elem_nodes_second[1];
    key_second.v3 = elem_nodes_second[2];
    key_second.v4 = elem_nodes_second[3];
 
    const bool same_key_first = (v.first.v1 == key_first.v1 &&
                                 v.first.v2 == key_first.v2 &&
                                 v.first.v3 == key_first.v3 &&
                                 v.first.v4 == key_first.v4);
    if (!same_key_first && quad_face.find(key_first) == quad_face.end()) {
      quad_inserts.push_back({key_first, v.second});
    }
    const bool same_key_second = (v.first.v1 == key_second.v1 &&
                                  v.first.v2 == key_second.v2 &&
                                  v.first.v3 == key_second.v3 &&
                                  v.first.v4 == key_second.v4);
    if (!same_key_second && quad_face.find(key_second) == quad_face.end()) {
      quad_inserts.push_back({key_second, v.second});
    }
  }
  for (const auto & entry : quad_inserts) quad_face.insert(entry);
}
 
template<class K> inline
void assign_counter_tags(hashmap::unordered_map<K, intersection_tags>& map, const MPI_Comm comm)
{
  size_t dsize = map.size();
  vector<K> key_vec(dsize);
  vector<intersection_tags> value_vec(dsize);
 
  size_t idx = 0;
  for (const auto& v : map) {
    key_vec[idx] = v.first;
    value_vec[idx] = v.second;
    idx++;
  }
 
  vector<int> perm, dest;
  emi_select_merge_destinations(key_vec, dsize, sizeof(K) + sizeof(intersection_tags), dest, comm,
                                "assign_counter_tags",
                                [](const K& key) { return hashmap::hash_ops<K>::hash(key); });
 
  commgraph<size_t> grph;
  grph.configure(dest, comm);
  size_t nrecv = sum(grph.rcnt);
 
  interval(perm, 0, dsize);
  binary_sort_copy(dest, perm);
 
  // fill send buffer and communicate
  vector<K> sbuff_key(dsize), rbuff_key(nrecv);
  vector<intersection_tags> sbuff_value(dsize), ibuff_value(dsize), rbuff_value(nrecv);
  for (size_t i = 0; i < dsize; i++) {
    sbuff_key[i] = key_vec[perm[i]];
    sbuff_value[i] = value_vec[perm[i]];
  }
  MPI_Exchange(grph, sbuff_key, rbuff_key, comm);
  MPI_Exchange(grph, sbuff_value, rbuff_value, comm);
 
  hashmap::unordered_map<K, intersection_tags> rmap;
  for (size_t i = 0; i < nrecv; i++) 
  {
    auto it = rmap.find(rbuff_key[i]);
    if (it != rmap.end()) 
    {
      intersection_tags& map_tags = it->second;
      intersection_tags& r_tags = rbuff_value[i];
 
      for (int r_tag : r_tags.tags) 
      {
        if (r_tag == -1) continue;
        bool found = false;
        for (int m_tag : map_tags.tags) 
        {
          if (m_tag == r_tag) {
            found = true;
            break;
          }
        }
        if (!found) 
        {
          for (size_t k = 0; k < MAX_INTERSECTIONS; ++k) {
            if (map_tags.tags[k] == -1) {
              map_tags.tags[k] = r_tag;
              break;
            }
          }
        }
      }
    }
    else
    {
      rmap.insert({ rbuff_key[i], rbuff_value[i] });
    }
  }
 
  for (size_t i = 0; i < nrecv; i++) {
    auto it = rmap.find(rbuff_key[i]);
    if (it != rmap.end()) rbuff_value[i] = it->second;
  }
 
  // sending the assigned data back to original rank
  grph.transpose();
  MPI_Exchange(grph, rbuff_value, ibuff_value, comm);
 
  for (size_t i = 0; i < dsize; i++) {
      auto it = map.find(sbuff_key[i]);
      if (it != map.end()) it->second = ibuff_value[i];
  }
}
 
template<class T, class S> inline
void compute_ptsdata_from_original_mesh(meshdata<T,S> & mesh, 
                                        const SF_nbr numbering,
                                        hashmap::unordered_map<T,T> & vertex2ptsdata, 
                                        hashmap::unordered_set<int> extra_tags,
                                        hashmap::unordered_set<int> intra_tags)
{
  MPI_Comm comm = mesh.comm;
  int      size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  const T* con = mesh.con.data();
  const SF::vector<T> & rnod = mesh.get_numbering(numbering);
 
  // Step 1: Gather all unique region tags associated with each vertex.
  // This map is local to each process and will store a list of tags for each vertex
  // based on the elements owned by the current process.
  hashmap::unordered_map<T, intersection_tags> map_index_to_tags;
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++)
  {
    T tag = mesh.tag[eidx];
    for (int n = mesh.dsp[eidx]; n < mesh.dsp[eidx+1];n++)
    {
      T gIndex = rnod[con[n]];
 
      auto it_new = map_index_to_tags.find(gIndex);
      if (it_new != map_index_to_tags.end())
      {
        // If vertex is already in the map, add the new tag if it's not already present.
        intersection_tags& tags = it_new->second;
        bool found = false;
        for(T t : tags.tags) {
          if(t == tag) {
            found = true;
            break;
          }
        }
        if(!found) {
          // This is an intracellular tag. Find an empty slot and add it.
          for(int i=0; i <MAX_INTERSECTIONS; ++i) {
            if (tags.tags[i] == -1) {
              tags.tags[i] = tag;
              break;
            }
          }
        }
      }
      else{
        // If vertex is not in the map, add it with the current tag.
        intersection_tags tags;
        tags.tags[0] = tag;
        map_index_to_tags.insert({gIndex,tags});
      }
    }
  }
 
  // Step 2: Globalize the tag information.
  // After this call, `map_index_to_tags` on every process will contain the
  // complete list of unique tags for every vertex in the entire mesh.
  assign_counter_tags(map_index_to_tags, comm);
 
  // Step 3: Classify each vertex based on its global list of tags.
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++)
  {
    for (int n = mesh.dsp[eidx]; n < mesh.dsp[eidx+1];n++)
    {
      T gIndex = rnod[con[n]];
      intersection_tags& vec_tags = map_index_to_tags[gIndex];
      int count_intra_tags = 0;
      int count_extra_tags = 0;
 
      // Count the number of unique intracellular and extracellular tags for the vertex.
      int count_tags = 0;
      for(T t : vec_tags.tags) {
        if(intra_tags.count(t)) count_intra_tags++;
        if(extra_tags.count(t)) count_extra_tags++;
      }
 
      // Apply classification rules based on the tag counts.
      if(count_extra_tags>=1 && count_intra_tags==0){
        vertex2ptsdata[gIndex] = 0; // Interior to an extracellular region
      }else if(count_extra_tags==0 && count_intra_tags==1){
        vertex2ptsdata[gIndex] = 0; // Interior to a myocyte
      }else if(count_extra_tags>=1 && count_intra_tags==1){
        vertex2ptsdata[gIndex] = 1; // On a membrane (1 myocyte, 1 extracellular)
      }else if(count_extra_tags==0 && count_intra_tags==2){
        vertex2ptsdata[gIndex] = 2; // On a gap junction (2 myocytes)
      }else if(count_extra_tags>=1 && count_intra_tags==2){
        vertex2ptsdata[gIndex] = 3; // On a complex junction (2 myocytes, 1 extracellular)
      }
      else if(count_intra_tags>2){
        // Error for unhandled cases.
        fprintf(stderr, "More than two intracellular are connected %lld. Tags:",
                static_cast<long long>(gIndex));
        for(T tag : vec_tags.tags) if(tag != -1) fprintf(stderr, " %lld", static_cast<long long>(tag));
        fprintf(stderr, "\n");
        EXIT(1);
      }
      else{
        // Error for unhandled cases.
        fprintf(stderr, "WARN: unhandled case in ptsData computation for gIndex %lld. Tags:",
                static_cast<long long>(gIndex));
        for(T tag : vec_tags.tags) if(tag != -1) fprintf(stderr, " %lld", static_cast<long long>(tag));
        fprintf(stderr, "\n");
        EXIT(1);
      }
    }
  }
}
 
 
template <class T, class S, class V, class emi_index_rank>
inline void assemble_map_both_to_unique(SF::abstract_matrix<T, S>&                                               op,
                                        const hashmap::unordered_map<mesh_int_t,
                                                                     std::pair<emi_index_rank, emi_index_rank>>& map,
                                        const SF::meshdata<mesh_int_t, V>&                                       unique_mesh,
                                        const SF::meshdata<mesh_int_t, V>&                                       both_mesh)
{
  // Recover global row/column ownership so locally assembled entries can be
  // inserted directly with global matrix indices.
  int rank;
  MPI_Comm_rank(both_mesh.comm, &rank);
 
  SF::vector<long int> layout_both;
  SF::layout_from_count<long int>(both_mesh.l_numelem, layout_both, both_mesh.comm);
 
  SF::vector<long int> layout_unique;
  SF::layout_from_count<long int>(unique_mesh.l_numelem, layout_unique, unique_mesh.comm);
 
  SF::vector<SF_int> row_idx(1), col_idx(1);
  SF::dmat<SF_real>  ebuff(1, 1);
 
  // Assemble one restriction row per locally owned unique face.
  for (const auto& [local_unique_idx, both_pair] : map) {
    const auto& first_both  = both_pair.first;
    const auto& second_both = both_pair.second;
 
    // The unique-face row is owned by this rank, so its global row index comes
    // from the unique-face layout plus the local unique-face element id.
    T global_unique_row = layout_unique[rank] + local_unique_idx;
    row_idx[0]          = global_unique_row;
 
    // Resolve the candidate both-face columns and discard invalid sides.
    bool first_valid  = (first_both.index >= 0 && first_both.rank >= 0);
    bool second_valid = (second_both.index >= 0 && second_both.rank >= 0);
 
    T global_both_col_first  = -1;
    T global_both_col_second = -1;
    if (first_valid) {
      global_both_col_first = layout_both[first_both.rank] + first_both.index;
    }
    if (second_valid) {
      global_both_col_second = layout_both[second_both.rank] + second_both.index;
    }
 
    // If both sides collapse to the same global both-face column, do not insert
    // the same contribution twice; treat it as a single-sided restriction row.
    if (first_valid && second_valid && global_both_col_first == global_both_col_second) {
      second_valid = false;
    }
 
    // Count how many distinct both-face columns contribute to this unique face.
    int count = 0;
    if (first_valid) count++;
    if (second_valid) count++;
 
    if (count == 0) continue;  // No data to average
 
    // Use arithmetic averaging when both sides exist; otherwise preserve the
    // single available value without halving it.
    double weight = 0.5;
    if (count == 1) weight = 1.0;
 
    ebuff.assign(1, 1, weight);
 
    // Insert the first contributing both-face column.
    if (first_valid) {
      col_idx[0] = global_both_col_first;
      op.set_values(row_idx, col_idx, ebuff.data(), false);
    }
 
    // Insert the second contributing both-face column when present.
    if (second_valid) {
      col_idx[0] = global_both_col_second;
      op.set_values(row_idx, col_idx, ebuff.data(), false);
    }
  }
 
  op.finish_assembly();
}
 
template<class T>
inline void restrict_to_membrane(vector<T> & v,
                                 hashmap::unordered_map<T,T> & dof2ptsData,
                                 const meshdata<mesh_int_t, mesh_real_t> & mesh)
{
  const SF::vector<SF_int> & rnod = mesh.get_numbering(SF::NBR_REF);
  hashmap::unordered_map<T,T> g2l;
  hashmap::unordered_map<T,T> l2g;
  for(size_t i=0; i<rnod.size(); i++){
    g2l[rnod[i]] = i;
    l2g[i] = rnod[i];
  }
 
  size_t widx = 0;
 
  for(size_t i=0; i<v.size(); i++){
    T g = l2g[v[i]];
 
    if (dof2ptsData[g] > 0) { // dof2ptsData requires global index
      v[widx++] = v[i];
    }
  }
 
  v.resize(widx);
}
 
}
 
#endif
#endif