emi.cc

Go to the documentation of this file.

//
// 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
 
#include "petsc_utils.h"
#include "physics_types.h"
#include "timers.h"
#include "stimulate.h"
#include "electric_integrators.h"
#include "SF_init.h" // for SF::init_xxx()
#include "emi.h"
#include <array>
 
#ifdef WITH_CALIPER
#include "caliper/cali.h"
#else
#include "caliper_hooks.h"
#endif
 
namespace opencarp {
 
void log_mesh_local_element_ranges(const sf_mesh& emi_mesh,
                                   const sf_mesh& emi_surfmesh_w_counter_face,
                                   const sf_mesh& emi_surfmesh_unique_face)
{
  int rank = 0;
  int comm_size = 0;
  MPI_Comm_rank(emi_surfmesh_w_counter_face.comm, &rank);
  MPI_Comm_size(emi_surfmesh_w_counter_face.comm, &comm_size);
 
  const unsigned long long local_emi_elems =
      static_cast<unsigned long long>(emi_mesh.l_numelem);
  const unsigned long long local_both_face_elems =
      static_cast<unsigned long long>(emi_surfmesh_w_counter_face.l_numelem);
  const unsigned long long local_unique_face_elems =
      static_cast<unsigned long long>(emi_surfmesh_unique_face.l_numelem);
 
  std::vector<unsigned long long> all_emi_elems;
  std::vector<unsigned long long> all_both_face_elems;
  std::vector<unsigned long long> all_unique_face_elems;
  if (rank == 0) {
    all_emi_elems.resize(comm_size, 0);
    all_both_face_elems.resize(comm_size, 0);
    all_unique_face_elems.resize(comm_size, 0);
  }
 
  MPI_Gather(&local_emi_elems, 1, MPI_UNSIGNED_LONG_LONG,
             rank == 0 ? all_emi_elems.data() : nullptr, 1, MPI_UNSIGNED_LONG_LONG,
             0, emi_surfmesh_w_counter_face.comm);
  MPI_Gather(&local_both_face_elems, 1, MPI_UNSIGNED_LONG_LONG,
             rank == 0 ? all_both_face_elems.data() : nullptr, 1, MPI_UNSIGNED_LONG_LONG,
             0, emi_surfmesh_w_counter_face.comm);
  MPI_Gather(&local_unique_face_elems, 1, MPI_UNSIGNED_LONG_LONG,
             rank == 0 ? all_unique_face_elems.data() : nullptr, 1, MPI_UNSIGNED_LONG_LONG,
             0, emi_surfmesh_w_counter_face.comm);
 
  if (rank != 0) return;
 
  const auto print_min_max = [](const char* label, const std::vector<unsigned long long>& counts) {
    if (counts.empty()) return;
 
    unsigned long long min_val = counts[0];
    unsigned long long max_val = counts[0];
    int min_rank = 0;
    int max_rank = 0;
 
    for (int r = 1; r < static_cast<int>(counts.size()); ++r) {
      if (counts[r] < min_val) {
        min_val = counts[r];
        min_rank = r;
      }
      if (counts[r] > max_val) {
        max_val = counts[r];
        max_rank = r;
      }
    }
 
    log_msg(NULL, 0, 0, "  %s: \n\t\t min=%llu on rank=%d, \n\t\t max=%llu on rank=%d\n",
            label, min_val, min_rank, max_val, max_rank);
  };
  log_msg(NULL, 0, 0, "\n**********************************");
  log_msg(NULL, 0, 0, "min/max number of local-element ranges:");
  print_min_max("emi_mesh", all_emi_elems);
  print_min_max("emi_surfmesh_w_counter_face", all_both_face_elems);
  print_min_max("emi_surfmesh_unique_face", all_unique_face_elems);
  log_msg(NULL, 0, 0, "**********************************");
}
 
void set_elec_tissue_properties_emi_volume(MaterialType* mtype, hashmap::unordered_set<int> &extra_tags, hashmap::unordered_set<int> &intra_tags, FILE_SPEC logger)
{
  CALI_CXX_MARK_FUNCTION;
  MaterialType *m = mtype;
  
  // initialize random conductivity fluctuation structure with PrM values
  m->regions.resize(param_globals::num_gregions);
 
  const char* grid_name = "emi_grid_domain";
  log_msg(logger, 0, 0, "Setting up %s tissue poperties for %d regions ..", grid_name,
          param_globals::num_gregions);
 
  char buf[64];
  RegionSpecs* reg = m->regions.data();
 
  // default tags for extra and intra domains
  hashmap::unordered_set<int> extra_tags_default = extra_tags;
  hashmap::unordered_set<int> intra_tags_default = intra_tags;
 
  for (size_t i=0; i<m->regions.size(); i++)
  {
    for (int j=0;j<param_globals::gregion[i].num_IDs;j++)
    {
      int tag = param_globals::gregion[i].ID[j];
 
      // removed tags explicitly defined by user in param_globals::gregion[i>1]
      if(extra_tags_default.find(tag) != extra_tags_default.end())
        extra_tags_default.erase(tag);
 
      if(intra_tags_default.find(tag) != intra_tags_default.end())
        intra_tags_default.erase(tag);
    }
  }
 
  for (size_t i=0; i<m->regions.size(); i++, reg++)  
  {
    if(!strcmp(param_globals::gregion[i].name, ""))  {
      snprintf(buf, sizeof buf, ", gregion_%d", int(i));
      param_globals::gregion[i].name = dupstr(buf);
    }
 
    // copy metadata into RegionSpecs
    reg->regname  = strdup(param_globals::gregion[i].name);
    reg->regID    = i;
 
    if(i==0) // default Extracellular region
      reg->nsubregs = extra_tags_default.size();
    if(i==1) // default intracellular region
      reg->nsubregs = intra_tags_default.size();
    if(i>1)  // optional: rest of other param_globals::gregion[i>1] defined by user 
      reg->nsubregs = param_globals::gregion[i].num_IDs;
 
    if(!reg->nsubregs)
      reg->subregtags = NULL;
    else  
    {
      reg->subregtags  = new int[reg->nsubregs];
 
      if(i==0){
        int j = 0;
        for (int tag : extra_tags_default) {
          reg->subregtags[j] = tag;
          j++;
        }
      }
      else if(i==1){
        int j = 0;
        for (int tag : intra_tags_default) {
          reg->subregtags[j] = tag;
          j++;
        }
      }
      else{
        for (int j=0;j<reg->nsubregs;j++)
          reg->subregtags[j] = param_globals::gregion[i].ID[j]; // explicit tags defined by user
      }
    }
 
    // describe material in given region
    elecMaterial *emat  = new elecMaterial();
    emat->material_type = ElecMat;
 
    // Isotropic conductivity is considered in EMI model.
    emat->InVal[0]   = param_globals::gregion[i].g_bath;
    emat->InVal[1]   = param_globals::gregion[i].g_bath;
    emat->InVal[2]   = param_globals::gregion[i].g_bath;
 
    emat->ExVal[0]   = param_globals::gregion[i].g_bath;
    emat->ExVal[1]   = param_globals::gregion[i].g_bath;
    emat->ExVal[2]   = param_globals::gregion[i].g_bath;
 
    emat->BathVal[0] = param_globals::gregion[i].g_bath;
    emat->BathVal[1] = param_globals::gregion[i].g_bath;
    emat->BathVal[2] = param_globals::gregion[i].g_bath;
 
    // convert units from S/m -> mS/um
    for (int j=0; j<3; j++) {
      emat->InVal[j]   *= 1e-3 * param_globals::gregion[i].g_mult;
      emat->ExVal[j]   *= 1e-3 * param_globals::gregion[i].g_mult;
      emat->BathVal[j] *= 1e-3 * param_globals::gregion[i].g_mult;
    }
    reg->material = emat;
  }
  
  if(strlen(param_globals::gi_scale_vec))
  {
    mesh_t mt = emi_msh;
    sf_mesh & mesh = get_mesh(mt);
    const char* file = param_globals::gi_scale_vec;
 
    size_t num_file_entries = SF::root_count_ascii_lines(file, mesh.comm);
 
    if(num_file_entries != mesh.g_numelem)
      log_msg(0,4,0, "%s warning: number of %s conductivity scaling entries does not match number of elements!",
              __func__, get_mesh_type_name(mt));
 
    // set up parallel element vector and read data
    sf_vec *escale;
    SF::init_vector(&escale, get_mesh(mt), 1, sf_vec::elemwise);
    escale->read_ascii(param_globals::gi_scale_vec);
 
    if(get_size() > 1) {
      // set up element vector permutation and permute
      if(get_permutation(mt, ELEM_PETSC_TO_CANONICAL, 1) == NULL) {
        register_permutation(mt, ELEM_PETSC_TO_CANONICAL, 1);
      }
      SF::scattering & sc = *get_permutation(mt, ELEM_PETSC_TO_CANONICAL, 1);
      sc(*escale, false);
    }
 
    // copy data into SF::vector
    SF_real* p = escale->ptr();
    m->el_scale.assign(p, p + escale->lsize());
    escale->release_ptr(p);
  }
}
 
void parabolic_solver_emi::init()
{
  CALI_CXX_MARK_FUNCTION;
  double t0, t1, dur;
  get_time(t0);
  stats.init_logger("par_stats.dat");
 
  // Create/initialise linear solver object: PETSc gets constructed with default settings
  SF::init_solver(&lin_solver);
 
  // EMI mesh: After decoupling the interfaces, we obtain a volumetric mesh with new dofs for all vertices.
  sf_mesh & emi_mesh = get_mesh(emi_msh);
  // Surface mesh called emi_surface_counter_msh:
  // This mesh contains only the interface faces, where the indices correspond
  // to the vertices of the original input mesh.
  // Each face has a corresponding counter-face, and they are identified by different tag numbers in each rank.
  sf_mesh & emi_surfmesh_w_counter_face = get_mesh(emi_surface_counter_msh);
  sf_mesh & emi_surfmesh_one_side = get_mesh(emi_surface_msh);
  sf_mesh & emi_surfmesh_unique_face = get_mesh(emi_surface_unique_face_msh);
 
  int max_row_entries_emi = max_nodal_edgecount(emi_mesh);
  int max_row_entries_emi_surfmesh = max_nodal_edgecount(emi_surfmesh_w_counter_face);
 
  int rank;
  MPI_Comm_rank(emi_surfmesh_w_counter_face.comm, &rank);
 
  sf_vec::ltype alg_type = sf_vec::algebraic;
  sf_vec::ltype alg_surface_type = sf_vec::elemwise;
 
  int dpn = 1;
 
  //-----------------------------------------------------------------
  // setup vectors
  //-----------------------------------------------------------------
  SF::init_vector(&ui,       emi_mesh, dpn, alg_type);
  SF::init_vector(&dui,      emi_mesh, dpn, alg_type);
  SF::init_vector(&ui_pre,   emi_mesh, dpn, alg_type);
  SF::init_vector(&Irhs,     emi_mesh, dpn, alg_type);
  SF::init_vector(&Iij_stim, emi_mesh, dpn, alg_type);
  SF::init_vector(&Iij_temp, emi_mesh, dpn, alg_type);
  SF::init_vector(&vb_both_face, emi_surfmesh_w_counter_face, dpn, alg_surface_type);
  SF::init_vector(&vb_unique_face, emi_surfmesh_unique_face, dpn, alg_surface_type);
  SF::init_vector(&Ib_both_face, emi_surfmesh_w_counter_face, dpn, alg_surface_type);
  SF::init_vector(&Ib_unique_face, emi_surfmesh_unique_face, dpn, alg_surface_type);
 
  //-----------------------------------------------------------------
  // initialize matrices (LHS, K, M_{emi mesh}, M_{surface mesh})
  //-----------------------------------------------------------------
  SF::init_matrix(&lhs_emi);
  SF::init_matrix(&stiffness_emi);
  SF::init_matrix(&mass_emi);
  SF::init_matrix(&mass_surf_emi);
 
  lhs_emi->init(emi_mesh, dpn, dpn, max_row_entries_emi*3);
  stiffness_emi->init(emi_mesh, dpn, dpn, max_row_entries_emi*2);
  mass_emi->init(emi_mesh, dpn, dpn, max_row_entries_emi*3);
  mass_surf_emi->init(emi_mesh, dpn, dpn, max_row_entries_emi*3);
  //---------------------------------------------------------------------
  // initialize operator matrices B, Bs, Vm_myocyte_emi
  //---------------------------------------------------------------------
  mesh_int_t M = emi_surfmesh_w_counter_face.g_numelem;
  mesh_int_t N = emi_mesh.g_numpts;
  mesh_int_t m = emi_surfmesh_w_counter_face.l_numelem;
  mesh_int_t m_one_side = emi_surfmesh_one_side.l_numelem;
  mesh_int_t M_one_side = emi_surfmesh_one_side.g_numelem; 
  mesh_int_t m_unique_face = emi_surfmesh_unique_face.l_numelem;
  mesh_int_t M_unique_face = emi_surfmesh_unique_face.g_numelem;
  mesh_int_t n = ui->lsize();
 
  log_msg(NULL, 0, 0, "\n**********************************");
  log_msg(NULL, 0, 0, "#elements of emi surfmesh unique face: %zu", emi_surfmesh_unique_face.g_numelem);
  log_msg(NULL, 0, 0, "#elements of emi surfmesh one side: %zu", emi_surfmesh_one_side.g_numelem);
  log_msg(NULL, 0, 0, "#elements of emi surfmesh: %zu", emi_surfmesh_w_counter_face.g_numelem);
  log_msg(NULL, 0, 0, "#elements of emi mesh: %zu", emi_mesh.g_numelem);
  log_msg(NULL, 0, 0, "#dofs for emi_mesh: %zu", emi_mesh.g_numpts);
  log_msg(NULL, 0, 0, "#max_row_entries_emi: %zu", max_row_entries_emi);
  log_msg(NULL, 0, 0, "**********************************\n");
  log_mesh_local_element_ranges(emi_mesh, emi_surfmesh_w_counter_face, emi_surfmesh_unique_face);
 
  SF::vector<long int> layout;
  SF::layout_from_count<long int>(emi_surfmesh_w_counter_face.l_numelem, layout, emi_surfmesh_w_counter_face.comm);
  mesh_int_t m_l = layout[rank];
  mesh_int_t n_l = emi_mesh.pl.algebraic_layout()[rank];
  const SF::vector<mesh_int_t> & alg_nod_surface = emi_surfmesh_w_counter_face.pl.algebraic_nodes();
 
  SF::vector<long int> layout_one_side;
  SF::layout_from_count<long int>(emi_surfmesh_one_side.l_numelem, layout_one_side, emi_surfmesh_one_side.comm);
  mesh_int_t m_one_side_l = layout_one_side[rank];
 
  SF::vector<long int> layout_unique_face;
  SF::layout_from_count<long int>(emi_surfmesh_unique_face.l_numelem, layout_unique_face, emi_surfmesh_unique_face.comm);
  mesh_int_t m_unique_face_l = layout_unique_face[rank];
 
  // To support a higher number of possible intersections between faces (e.g., multiple contacts on the membrane),
  // such as when one myocyte is surrounded by more than one extracellular region,
  // we slightly overallocate the maximum number of row entries (max_row_entries)
  // for the matrices B, Bi, BsM, and Vm_myocyte_emi to ensure sufficient storage for all coupling cases.
  SF::init_matrix(&B);
  B->init(M, N, m, n, m_l, max_row_entries_emi*2);
  B->zero();
  SF::init_matrix(&Bi);
  Bi->init(M, N, m, n, m_l, max_row_entries_emi*2);
  Bi->zero();
  SF::init_matrix(&BsM);
  BsM->init(N, M, n, m, n_l, max_row_entries_emi*2);
  BsM->zero();
  SF::init_matrix(&Vm_myocyte_emi); 
  Vm_myocyte_emi->init(N, N, n, n, n_l, max_row_entries_emi*2);
  Vm_myocyte_emi->zero();
 
  // Build the direct unique-face <-> both-face transfer operators
  SF::construct_direct_unique_both_operators(operator_unique_to_both_faces,
                                             operator_both_to_unique_face,
                                             map_elem_uniqueFace_to_elem_bothface,
                                             map_elem_uniqueFace_to_elem_oneface,
                                             vec_both_to_one_face,
                                             emi_surfmesh_w_counter_face,
                                             emi_surfmesh_unique_face,
                                             max_row_entries_emi,
                                             dpn,
                                             alg_surface_type);
 
  // initialize Interpolation operator by B and Scatter operator by Bs
  SF::assemble_restrict_operator(*B, *Bi, *BsM, elemTag_surface_w_counter_mesh, map_vertex_tag_to_dof_petsc, line_face, tri_face, quad_face, emi_surfmesh_w_counter_face, emi_mesh, UM2_to_CM2);
 
  // DEBUG: Check mesh sizes and mappings
  #ifdef EMI_DEBUG_MESH
  {
    int local_rank;
    MPI_Comm_rank(emi_mesh.comm, &local_rank);
    
    fprintf(stderr, "RANK %d MESH SIZES: one_side=%zu, counter=%zu, unique=%zu\n",
            local_rank, emi_surfmesh_one_side.l_numelem, 
            emi_surfmesh_w_counter_face.l_numelem, emi_surfmesh_unique_face.l_numelem);
    
    fprintf(stderr, "RANK %d MAP SIZES: uniqueFace_to_oneface=%zu\n",
            local_rank, map_elem_uniqueFace_to_elem_oneface.size());
    fflush(stderr);
  }
  #endif
 
  //-----------------------------------------------------------------
  // setup the ionic current
  //-----------------------------------------------------------------
  // get initial value of the vm and Iion on the face from the ionicOnFace class
  sf_vec* vb_ptr   = get_data(vm_emi_itf);
  sf_vec* Ib_ptr = get_data(iion_emi_itf);
 
  if(!(vb_ptr != NULL && Ib_ptr != NULL)) {
    log_msg(0,5,0, "%s error: global Vb and Ib vectors not properly set up! Ionics seem invalid! Aborting!",
            __func__);
    EXIT(1);
  }
 
  SF::init_vector(&vb, vb_ptr);
  vb->shallow_copy(*vb_ptr);
  SF::init_vector(&Ib, Ib_ptr);
  Ib->shallow_copy(*Ib_ptr);
 
  parab_tech = static_cast<parabolic_solver_emi::parabolic_t>(param_globals::parab_solve_emi);
 
  dur = timing(t1, t0);
}
 
void parabolic_solver_emi::rebuild_matrices(MaterialType* mtype, limpet::MULTI_IF & miif, SF::vector<stimulus> & stimuli, FILE_SPEC logger)
{
  CALI_CXX_MARK_FUNCTION;
  double start, end, period;
  get_time(start);
  double t0, t1, dur;
  mass_integrator mass_integ;
  mass_integrator mass_integ_emi;
  int dpn = 1;
 
  int            log_flag = param_globals::output_level > 1 ? ECHO : 0;
  MaterialType & mt       = mtype[0];
  const bool have_dbc = have_dbc_stims(stimuli);
 
  double Dt = user_globals::tm_manager->time_step;
 
  cond_t    condType = intra_cond;
  sf_mesh & mesh     = get_mesh(emi_msh);
  sf_mesh & emi_surfmesh_w_counter_face    = get_mesh(emi_surface_counter_msh);
 
  // set material and conductivity type
  set_cond_type(mt, condType);
 
  // assemble EMI matrices
  // fill the EMI system
  {
    // emi-> stiffness
    log_msg(NULL, 0, 0, "assemble stiffness matrix");
    get_time(t0);
    elec_stiffness_integrator stfn_integ_emi(mt);
    stiffness_emi->zero();
    SF::assemble_matrix(*stiffness_emi, mesh, stfn_integ_emi);
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "Computed parabolic stiffness matrix in %.5f seconds.", float(dur));
 
    log_msg(NULL, 0, 0, "assemble mass matrix on the volumetric mesh");
    get_time(t0);
    mass_integrator mass_integ;
    mass_emi->zero();
    SF::assemble_matrix(*mass_emi, mesh, mass_integ);
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "Computed volumetric mass matrix in %.5f seconds.", float(dur));
 
    log_msg(NULL, 0, 0, "assemble LHS matrix and mass matrix on the surface mesh");
    get_time(t0);
    lhs_emi->zero();
    // lhs_emi = (CmM/dt + K)
    SF::assemble_lhs_emi(*lhs_emi, *mass_surf_emi, mesh, emi_surfmesh_w_counter_face, map_vertex_tag_to_dof_petsc, line_face, tri_face, quad_face, stfn_integ_emi, mass_integ_emi, -1., UM2_to_CM2 / Dt);
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "Computed parabolic mass matrix in %.5f seconds.", float(dur));
  }
 
 
  bool same_nonzero = false;
 
  // set boundary conditions
  if(have_dbc) {
    log_msg(logger,0,log_flag, "lhs matrix enforcing Dirichlet boundaries.");
    get_time(t0);
 
    if(dbc == nullptr)
      dbc = new dbc_manager(*lhs_emi, stimuli);
    else
      dbc->recompute_dbcs();
 
    dbc->enforce_dbc_lhs();
 
    dur = timing(t1,t0);
    log_msg(logger,0,log_flag, "lhs matrix Dirichlet enforcing done in %.5f seconds.", float(dur));
  }
  else {
    log_msg(logger,0,ECHO, "without enforcing Dirichlet boundaries on the lhs matrix!");
    // we are dealing with a singular system
    phie_mat_has_nullspace = true;
  }
 
  set_dir(INPUT);
  get_time(t0);
  setup_linear_solver(logger);
  dur = timing(t1,t0);
  log_msg(logger,0,log_flag, "Initializing parabolic solver in %.5f seconds.", float(dur));
  set_dir(OUTPUT);
 
  period = timing(end, start);
}
 
void parabolic_solver_emi::setup_linear_solver(FILE_SPEC logger)
{
  tol    = param_globals::cg_tol_parab;
  max_it = param_globals::cg_maxit_parab;
 
  std::string default_opts;
  std::string solver_file;
  solver_file = param_globals::parab_options_file;
  if (param_globals::flavor == std::string("ginkgo")) {
    default_opts = std::string(
        R"(
{
    "type": "solver::Cg",
    "preconditioner": {
        "type": "solver::Multigrid",
        "min_coarse_rows": 8,
      "max_levels": 16,
        "default_initial_guess": "zero",
        "mg_level": [
            {
                "type": "multigrid::Pgm",
                "deterministic": false
            }
        ],
        "coarsest_solver": {
             "type": "preconditioner::Schwarz",
             "local_solver": {
                 "type": "preconditioner::Jacobi"
             }
        },
        "criteria": [
            {
                "type": "Iteration",
                "max_iters": 1
            }
        ]
    },
    "criteria": [
        {
            "type": "Iteration",
            "max_iters": 100
        },
        {
            "type": "ResidualNorm",
            "reduction_factor": 1e-4
        }
    ]
}
      )");
  } else if (param_globals::flavor == std::string("petsc")) {
    default_opts = std::string("-ksp_type cg -pc_type gamg -options_left");
  }
  lin_solver->setup_solver(*lhs_emi, tol, max_it * 100, param_globals::cg_norm_parab,
                           "parabolic PDE", phie_mat_has_nullspace, logger, solver_file.c_str(),
                           default_opts.c_str());
}
 
void parabolic_solver_emi::solve()
{
  switch (parab_tech) {
    case SEMI_IMPLICIT:  solve_semiImplicit();   break;
  }
}
 
void parabolic_solver_emi::solve_semiImplicit()
{
  double t0,t1;
  get_time(t0);
 
  if(dbc != nullptr){  
    CALI_MARK_BEGIN("apply_dbc_rhs");
    dbc->enforce_dbc_rhs(*ui);
    CALI_MARK_END("apply_dbc_rhs");
  }
 
  *ui_pre = *ui;
 
  // -K*u (K is assembled as -K)
  CALI_MARK_BEGIN("stiff_mat");
  stiffness_emi->mult(*ui, *Iij_temp);
  // rhs = Iij + K*u
  *Irhs -= *Iij_temp;
  CALI_MARK_END("stiff_mat");
 
  // add volumetric stimulus currents (I_ex, I_in)
  // rhs = Iij + K*u + M*Iij_stim
  if (Iij_stim->mag() > 0.0) {
    CALI_MARK_BEGIN("stim_application");
    mass_emi->mult(*Iij_stim, *Iij_temp);
    *Irhs -= *Iij_temp;
    CALI_MARK_BEGIN("stim_application");
  }
 
  // rhs = -(Iij + K*u + M*Iij_stim)
  CALI_MARK_BEGIN("rhs_update");
  (*Irhs) *= (-1.0);
  CALI_MARK_END("rhs_update");
 
  // compute step
  CALI_MARK_BEGIN("linear_solve");
  (*lin_solver)(*dui, *Irhs);    
  CALI_MARK_END("linear_solve");
 
  // logfile for solver
  if(lin_solver->reason < 0) {
    log_msg(0, 5, 0,"%s solver diverged. Reason: %s.", lin_solver->name.c_str(),
            petsc_get_converged_reason_str(lin_solver->reason));
    EXIT(1);
  }
 
  // update solution:: ui_pre = ui_pre + dui
  CALI_MARK_BEGIN("sol_update");
  ui_pre->add_scaled(*dui, 1.0);
  *ui *=0;
  ui->add_scaled(*ui_pre, 1.0);
  CALI_MARK_END("sol_update");
 
  // We need to enforce DBCs again after the solution vector was updated.
  // Otherwise, matrix/solver tolerances allow small nonzero values in dui to accumulate in ui.
  if(dbc != nullptr){  
    CALI_MARK_BEGIN("apply_dbc_rhs");
    dbc->enforce_dbc_rhs(*ui);
    CALI_MARK_END("apply_dbc_rhs");
  }
 
  // treat solver statistics
  stats.slvtime += timing(t1, t0);
  stats.update_iter(lin_solver->niter);
}
 
void tags_onFace( hashmap::unordered_map<SF::tuple<mesh_int_t>, 
                      std::pair<SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>,
                      SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>>> & line_face,
                      hashmap::unordered_map<SF::triple<mesh_int_t>, 
                      std::pair<SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>,
                      SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>>> & tri_face,
                      hashmap::unordered_map<SF::quadruple<mesh_int_t>, 
                      std::pair<SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>, 
                      SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>>> & quad_face,
                      hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, mesh_int_t> & map_vertex_tag_to_dof,
                      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,
                      std::vector<std::string> & tags_data)
{
  sf_mesh & emi_surfmesh_w_counter_face = get_mesh(emi_surface_counter_msh);
 
  const SF::vector<mesh_int_t> & tags = emi_surfmesh_w_counter_face.tag;
 
  const SF::vector<mesh_int_t> & rnod = emi_surfmesh_w_counter_face.get_numbering(SF::NBR_REF);
  hashmap::unordered_map<mesh_int_t,mesh_int_t> l2g;
 
  for(size_t i=0; i<rnod.size(); i++){
    l2g[i] = rnod[i];
  }
 
  for(size_t eidx=0; eidx<emi_surfmesh_w_counter_face.l_numelem; eidx++)
  {
    std::vector<int> elem_nodes;
    mesh_int_t tag = emi_surfmesh_w_counter_face.tag[eidx];
    for (int n = emi_surfmesh_w_counter_face.dsp[eidx]; n < emi_surfmesh_w_counter_face.dsp[eidx+1];n++)
    {
      mesh_int_t l_idx = emi_surfmesh_w_counter_face.con[n];
 
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(l2g[l_idx],tags[eidx]);
      mesh_int_t dof = map_vertex_tag_to_dof[Index_tag_old];
      elem_nodes.push_back(dof);
    }
 
    mesh_int_t tag_first = 0;
    mesh_int_t tag_second = 0;
    std::string result_first;
    std::string result_second;
    std::sort(elem_nodes.begin(),elem_nodes.end()); // make the node tuple order-independent
 
    // Extract all surface face pairs separating regions with different material or boundary tags (ionicFaces)
    if(elem_nodes.size()==2){
      SF::tuple<mesh_int_t> key;
 
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      std::pair<SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>,
                SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>> value = line_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
      result_first = std::to_string(tag_first) + ":" + std::to_string(tag_second);
      result_second = std::to_string(tag_second) + ":" + std::to_string(tag_first); 
    }
    else if(elem_nodes.size()==3){
      SF::triple<mesh_int_t> key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      std::pair<SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>,
                SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>>  value = tri_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
      result_first = std::to_string(tag_first) + ":" + std::to_string(tag_second);
      result_second = std::to_string(tag_second) + ":" + std::to_string(tag_first);
    }
    else if(elem_nodes.size()==4){
      SF::quadruple<mesh_int_t> key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      key.v4 = elem_nodes[3];
      std::pair<SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>,
                SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>> value = quad_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
      result_first = std::to_string(tag_first) + ":" + std::to_string(tag_second);
      result_second = std::to_string(tag_second) + ":" + std::to_string(tag_first);
    }
    // check if current element tag has a custom region ID
    tags_data.push_back(result_first);
    tags_data.push_back(result_second);
  }
}
 
void EMI::initialize()
{
  double t1, t2;
  get_time(t1);
 
  set_dir(OUTPUT);
 
  // open logger
  logger = f_open("emi.log", param_globals::experiment != 4 ? "w" : "r");
  const int verb = param_globals::output_level;
  // Mesh processing step: convert the input mesh into 
  // - an EMI mesh (with discontinuities at gap junctions and membranes) 
  // - and generate the corresponding surface mesh
  CALI_MARK_BEGIN("mesh_setup");
  setup_EMI_mesh();
 
  // setup mappings between extra and intra grids, algebraic and nodal,
  // and between PETSc and canonical orderings
  setup_mappings();
 
  // the ionicOnFace physics is currently triggered from inside the emi to have tighter
  // control over it
  ion.logger = logger;
 
  ion.set_surface_mesh_data(parab_solver.line_face,
                            parab_solver.tri_face,
                            parab_solver.quad_face,
                            parab_solver.map_vertex_tag_to_dof);
 
  // builds per-face adjacency tags
  std::vector<std::string> tags_data;                    
  tags_onFace(parab_solver.line_face,
              parab_solver.tri_face,
              parab_solver.quad_face,
              parab_solver.map_vertex_tag_to_dof,
              parab_solver.map_vertex_tag_to_dof_petsc,
              tags_data);
  CALI_MARK_END("mesh_setup");
 
  ion.set_tags_onFace(tags_data);
 
  ion.set_face_region_data(parab_solver.intra_tags, parab_solver.map_elem_uniqueFace_to_tags);
  ion.initialize();
 
  // set up tissue properties on the extra and intracellular domains
  set_elec_tissue_properties_emi_volume(mtype_vol, parab_solver.extra_tags, parab_solver.intra_tags, logger);
  region_mask(emi_msh, mtype_vol[0].regions, mtype_vol[0].regionIDs, true, "gregion_vol");
  
  // add electrics timer for time stepping, add to time stepper tool (TS)
  double global_time = user_globals::tm_manager->time;
  timer_idx = user_globals::tm_manager->add_eq_timer(global_time, param_globals::tend, 0,
                param_globals::dt, 0, "elec::ref_dt", "TS");
 
  // EMI stimuli setup
  CALI_MARK_BEGIN("stimulus_setup");
  param_globals::operator_splitting = 0; // as long as we don't use operator splitting, we default to 0 so we can keep the stimuli scaling mint
  setup_stimuli();
  CALI_MARK_END("stimulus_setup");
 
  // set up the linear equation systems. this needs to happen after the stimuli have been
  // set up, since we need boundary condition info
  CALI_MARK_BEGIN("solver_setup");
  setup_solvers();
  CALI_MARK_END("solver_setup");
  
  // balance electrodes, we may need the extracellular mass matrix
  balance_electrodes();
  // total current scaling
  scale_total_stimulus_current(stimuli, *parab_solver.mass_emi, *parab_solver.mass_surf_emi, logger);
 
  sf_mesh & emi_mesh = get_mesh(emi_msh);
  sf_mesh & emi_surfmesh_w_counter_face = get_mesh(emi_surface_counter_msh);
 
  // - Extract the resting potential from v_b, which is defined only on the membrane faces. 
  //   (we do not care about gap junction in this case)
  // - Project this resting potential from the surface mesh onto the EMI mesh.
  // - assign the corresponding resting potential on the memebrane to all DOFs 
  //   myocyte that is associated with those membrane faces. 
  //   (dofs on the interfaces of myocytes and all inner dofs located inside of myocytes).
  // - any DOF with tag number belongs to the extra_tags set, is assigned to zero.
  // - We provide a mapping between the vertices of the surface mesh(vertices of the original input mesh) 
  //   and those of the EMI mesh (with updated indices).
  //   Based on the face type (line_face, tri_face, or quad_face) and the tag numbers that 
  //   identify paired faces (face and counter-face), 
  //   we can determine whether a given face belongs to a membrane or a gap-junction region.
  // - Direct unique -> both mapping
  parab_solver.operator_unique_to_both_faces->mult(*parab_solver.vb, *parab_solver.vb_both_face);
  SF::assign_resting_potential_from_ionic_models_on_myocyte(*parab_solver.Vm_myocyte_emi,
                    parab_solver.vb_both_face,
                    parab_solver.elemTag_emi_mesh,
                    parab_solver.map_vertex_tag_to_dof_petsc, 
                    parab_solver.line_face, parab_solver.tri_face, parab_solver.quad_face, 
                    emi_surfmesh_w_counter_face, emi_mesh);  
 
  parab_solver.Vm_myocyte_emi->get_diagonal(*parab_solver.ui); // initialize all myocytes to the resting potential based on ionic models
  *parab_solver.vb_unique_face = *parab_solver.vb;
 
  CALI_MARK_BEGIN("output_setup");
  // prepare the electrics output. we skip it if we do post-processing
  if(param_globals::experiment != EXP_POSTPROCESS)
    setup_output();
  CALI_MARK_END("output_setup");
 
  this->initialize_time += timing(t2, t1);
}
 
void EMI::setup_mappings()
{
  bool emi_exits = mesh_is_registered(emi_msh);
  assert(emi_exits);
  const int dpn = 1;
 
  if(get_scattering(emi_msh, ALG_TO_NODAL, dpn) == NULL) 
  {
    log_msg(logger, 0, 0, "%s: Setting up intracellular algebraic-to-nodal scattering.", __func__);
    register_scattering(emi_msh, ALG_TO_NODAL, dpn);
  }
  if(get_permutation(emi_msh, PETSC_TO_CANONICAL, dpn) == NULL) 
  {
    log_msg(logger, 0, 0, "%s: Setting up intracellular PETSc to canonical permutation.", __func__);
    register_permutation(emi_msh, PETSC_TO_CANONICAL, dpn);
  }
}
 
void EMI::compute_step()
{
  double t1, t2;
  get_time(t1);
 
  // activation checking
  const double time      = user_globals::tm_manager->time,
               time_step = user_globals::tm_manager->time_step;
 
  const int verb = param_globals::output_level;
  // We are treating stimuli by type:
  //    - Potential stimuli (Phi_ex, GND_ex, Phi_ex_ol) are managed by a dbc_manager and are applied 
  //      to the left and right hand side of the equation.
  //    - Current stimuli (I_ex, I_in) are applied to the right hand side,
  //      while (I_tm) is applied to the vector Ib directly.
  CALI_MARK_BEGIN("apply_dbc_lhs");
  apply_dbc_stimulus();
  CALI_MARK_END("apply_dbc_lhs");
 
  // compute ionics update
  CALI_MARK_BEGIN("ion_compute");
  ion.compute_step();
  CALI_MARK_END("ion_compute");
 
  CALI_MARK_BEGIN("apply_stim");
  apply_current_stimulus();
  CALI_MARK_END("apply_stim");
 
  // convert Ib -> Irhs
  parab_solver.operator_unique_to_both_faces->mult(*parab_solver.Ib, *parab_solver.Ib_both_face);
  parab_solver.BsM->mult(*parab_solver.Ib_both_face, *parab_solver.Irhs);
  
  // solver parabolic system
  CALI_MARK_BEGIN("parab_solve");
  parab_solver.solve();
  {
    // v_b = B_i * u
    parab_solver.B->mult(*parab_solver.ui, *parab_solver.vb_both_face);
    // Direct both -> unique mapping
    parab_solver.operator_both_to_unique_face->mult(*parab_solver.vb_both_face, *parab_solver.vb_unique_face);
    *parab_solver.vb = *parab_solver.vb_unique_face;
  }
  CALI_MARK_END("parab_solve");
 
  if(user_globals::tm_manager->trigger(iotm_console)) {
    // output lin solver stats
    parab_solver.stats.log_stats(user_globals::tm_manager->time, false);
  }
  this->compute_time += timing(t2, t1);
 
  // since the traces have their own timing, we check for trace dumps in the compute step loop
  if(user_globals::tm_manager->trigger(iotm_trace))
    dump_trace( ion.miif, user_globals::tm_manager->time );
}
 
void EMI::output_step()
{
  double t1, t2;
  get_time(t1);
  
  output_manager.write_data();
 
  double curtime = timing(t2, t1);
  this->output_time += curtime;
 
  IO_stats.calls++;
  IO_stats.tot_time += curtime;
 
  if(user_globals::tm_manager->trigger(iotm_console))
    IO_stats.log_stats(user_globals::tm_manager->time, false);
}
 
void EMI::destroy()
{
  CALI_CXX_MARK_FUNCTION;
  // close logger
  f_close(logger);
 
  // close output files
  output_manager.close_files_and_cleanup();
 
  // destroy ionics
  ion.destroy();
}
 
void EMI::setup_stimuli()
{
  // initialize basic stim info data (used units, supported types, etc)
  init_stim_info();
  
  stimuli.resize(param_globals::num_stim);
  for (int i = 0; i < param_globals::num_stim; i++) {
    // construct new stimulus
    stimulus & s = stimuli[i];
 
    if (user_globals::using_legacy_stimuli)
      s.translate(i);
 
    // we associate to the EMI mesh. this is needed for the stim_phys and stim_electrode setups.
    s.associated_intra_mesh = emi_msh, s.associated_extra_mesh = emi_msh;
 
    s.setup(i);
 
    if (s.phys.type == Illum) {
      log_msg(0, MAX_LOG_LEVEL, ECHO, "Stimulus of type Illum (=6) is not implemented in EMI. Abort.");
      EXIT(EXIT_FAILURE);
    }
 
    // Depending on the stimulus type, we make sure to only stimulate the correct regions of the mesh:
    // Extracellular stimuli restrict to all DOFs of the extracellular region (including the extracellular side of the membrane).
    // Equivalent for intracellular stimuli.
    // Stimuli that act directly on the membrane restrict to DOFs with ptsData > 0, i.e. membrane, gap junctions, and complex junctions.
    if (is_extra(s.phys.type)) {
      SF::vector<mesh_int_t> extra_vertices;
      const sf_mesh&         mesh = get_mesh(s.associated_extra_mesh);
 
      // Gather vertices from emi_msh using extra_tags
      indices_from_region_tags(extra_vertices, mesh, parab_solver.extra_tags);
 
      // Restrict electrode vertices to extra region
      restrict_to_set(s.electrode.vertices, extra_vertices);
    } else if (s.phys.type == I_tm) {
      const sf_mesh&         mesh = get_mesh(emi_msh);
      SF::restrict_to_membrane(s.electrode.vertices, dof2ptsData, mesh);
    } else {
      SF::vector<mesh_int_t> intra_vertices;
      const sf_mesh&         mesh = get_mesh(s.associated_intra_mesh);
 
      // Gather vertices from emi_msh using intra_tags
      indices_from_region_tags(intra_vertices, mesh, parab_solver.intra_tags);
 
      // Restrict electrode vertices to intra region
      restrict_to_set(s.electrode.vertices, intra_vertices);
    }
 
    if (s.electrode.dump_vtx) {
      set_dir(OUTPUT);
      s.dump_vtx_file(i);
    }
 
    if(param_globals::stim[i].pulse.dumpTrace && get_rank() == 0) {
      set_dir(OUTPUT);
      s.pulse.wave.write_trace(s.name+".trc");
    }
 
  }
}
 
void EMI::apply_dbc_stimulus()
{
  parabolic_solver_emi& ps = parab_solver;
 
  // we check if the DBC layout changed, if so we recompute the matrix and the dbc_manager
  bool dbcs_have_updated = ps.dbc != nullptr && ps.dbc->dbc_update();
  bool time_not_final    = user_globals::tm_manager->d_time != user_globals::tm_manager->d_end;
 
  if (dbcs_have_updated && time_not_final) {
    parab_solver.rebuild_matrices(mtype_vol, *ion.miif, stimuli, logger);
  }
}
 
void EMI::apply_current_stimulus()
{
  parabolic_solver_emi& ps = parab_solver;
  ps.Iij_stim->set(0.0);
 
  // iterate over stimuli
  for(stimulus & s : stimuli) {
    if(s.is_active()) {
      switch (s.phys.type) {
        case I_tm: {
          apply_stim_to_vector(s, *ps.Iij_temp, true);
          ps.Bi->mult(*ps.Iij_temp, *ps.Ib_both_face);
          ps.operator_both_to_unique_face->mult(*ps.Ib_both_face, *ps.Ib_unique_face);
          ps.Ib->add_scaled(*ps.Ib_unique_face, -0.5); // scale by 0.5 to revert 2x scaling in Bi*Iij
        } break;
 
        case I_ex:
        case I_in: {
          apply_stim_to_vector(s, *ps.Iij_stim, true);
        } break;
 
        default: break;
      }
    }
  }
}
 
void EMI::balance_electrodes()
{
  for (int i = 0; i < param_globals::num_stim; i++) {
    if (param_globals::stim[i].crct.balance != -1) {
      int from = param_globals::stim[i].crct.balance;
      int to   = i;
 
      log_msg(NULL, 0, 0, "Balancing stimulus %d with %d %s-wise.", from, to,
              is_current(stimuli[from].phys.type) ? "current" : "voltage");
 
      stimulus& s_from = stimuli[from];
      stimulus& s_to   = stimuli[to];
 
      s_to.pulse = s_from.pulse;
      s_to.ptcl  = s_from.ptcl;
      s_to.phys  = s_from.phys;
      s_to.pulse.strength *= -1.0;
 
      if (s_from.phys.type == I_ex || s_from.phys.type == I_in) {
        // if from is total current, skip volume based adjustment of strength
        // otherwise, calling constant_total_stimulus_current() will undo the balanced scaling of to.pulse.strength
        // constant_total_stimulus_current() will do the scaling based on the volume
        if (!s_from.phys.total_current) {
          sf_mat& mass = *parab_solver.mass_emi;
          SF_real vol0 = get_volume_from_nodes(mass, s_from.electrode.vertices);
          SF_real vol1 = get_volume_from_nodes(mass, s_to.electrode.vertices);
 
          s_to.pulse.strength *= fabs(vol0 / vol1);
        }
      }
    }
  }
}
 
void EMI::scale_total_stimulus_current(SF::vector<stimulus>& stimuli,
                                       sf_mat&               mass_vol,
                                       sf_mat&               mass_surf,
                                       FILE_SPEC             logger)
{
  for (stimulus & s : stimuli){
    if(is_current(s.phys.type) && s.phys.total_current){
      switch (s.phys.type) {
        case I_in:
        case I_ex: {
          // compute affected volume in um^3
          SF_real vol = get_volume_from_nodes(mass_vol, s.electrode.vertices);
          // s->strength holds the total current in uA, compute current density in uA/cm^3
          // Theoretically, we don't need to scale the volume to cm^3 here since we later
          // multiply with the mass matrix and we get um^3 * uA/um^3 = uA.
          // However, for I_ex/I_in there is an additional um^3 to cm^3 scaling in phys.scale,
          // since I_ex/I_in is expected to be in uA/cm^3. Therefore, we need to compensate for that to arrive at uA later.
          assert(vol > 0);
          float scale = 1.e12 / vol;
 
          s.pulse.strength *= scale;
 
          log_msg(logger, 0, ECHO,
                  "%s [Stimulus %d]: current density scaled to %.4g uA/cm^3\n",
                  s.name.c_str(), s.idx, s.pulse.strength);
        } break;
 
        case I_tm: {
          // In the EMI model, I_tm only affects the membrane. Therefore, we compute
          // the affected membrane surface in um^2 using the membrane mass matrix.
          // The electrode vertices are resticted to the membrane during setup, hence this function returns a surface already.
          SF_real surf = get_volume_from_nodes(mass_surf, s.electrode.vertices);
 
          // convert to cm^2
          assert(surf > 0);
          surf /= 1.e8;
 
          // scale surface density now to result in correct total current
          s.pulse.strength /= surf;
          log_msg(logger, 0, ECHO,
                  "%s [Stimulus %d]: current density scaled to %.4g uA/cm^2\n",
                  s.name.c_str(), s.idx, s.pulse.strength);
        } break;
 
        default: break;
      }
    }
  }
}
 
// Assign a deterministic global element numbering for a surface mesh based on
// element type, tag, and sorted node ids. This only affects output ordering.
static void assign_deterministic_elem_numbering(sf_mesh & mesh)
{
  const int KEY_SIZE = 6; // type, tag, n1, n2, n3, n4
  int rank = 0, size = 0;
  MPI_Comm_rank(mesh.comm, &rank);
  MPI_Comm_size(mesh.comm, &size);
 
  auto make_key = [&](size_t i) {
    std::array<long long, KEY_SIZE> k;
    k.fill(-1);
    k[0] = (long long)mesh.type[i];
    k[1] = (long long)mesh.tag[i];
 
    int nn = 0;
    if (mesh.type[i] == SF::Line) nn = 2;
    else if (mesh.type[i] == SF::Tri) nn = 3;
    else if (mesh.type[i] == SF::Quad) nn = 4;
    else nn = 0;
 
    std::vector<long long> nodes;
    nodes.reserve(nn);
    size_t off = mesh.dsp[i];
    for (int j = 0; j < nn; j++) {
      nodes.push_back((long long)mesh.con[off + j]);
    }
    std::sort(nodes.begin(), nodes.end());
    for (int j = 0; j < (int)nodes.size(); j++) {
      k[2 + j] = nodes[j];
    }
    return k;
  };
 
  // Pack local keys
  std::vector<long long> local_keys(mesh.l_numelem * KEY_SIZE, -1);
  for (size_t i = 0; i < mesh.l_numelem; i++) {
    auto k = make_key(i);
    for (int j = 0; j < KEY_SIZE; j++) local_keys[i * KEY_SIZE + j] = k[j];
  }
 
  // Gather sizes
  std::vector<int> counts(size, 0), displs(size, 0);
  int local_count = (int)local_keys.size();
  MPI_Allgather(&local_count, 1, MPI_INT, counts.data(), 1, MPI_INT, mesh.comm);
  int total = 0;
  for (int r = 0; r < size; r++) {
    displs[r] = total;
    total += counts[r];
  }
 
  std::vector<long long> all_keys;
  if (rank == 0) all_keys.resize(total, -1);
  MPI_Gatherv(local_keys.data(), local_count, MPI_LONG_LONG,
              rank == 0 ? all_keys.data() : nullptr, counts.data(), displs.data(), MPI_LONG_LONG,
              0, mesh.comm);
 
  // Build sorted global keys on rank 0
  std::vector<std::array<long long, KEY_SIZE>> sorted_keys;
  if (rank == 0) {
    const int nkeys = total / KEY_SIZE;
    sorted_keys.resize(nkeys);
    for (int i = 0; i < nkeys; i++) {
      std::array<long long, KEY_SIZE> k;
      for (int j = 0; j < KEY_SIZE; j++) k[j] = all_keys[i * KEY_SIZE + j];
      sorted_keys[i] = k;
    }
    std::sort(sorted_keys.begin(), sorted_keys.end());
  }
 
  // Broadcast sorted keys
  int nkeys = 0;
  if (rank == 0) nkeys = (int)sorted_keys.size();
  MPI_Bcast(&nkeys, 1, MPI_INT, 0, mesh.comm);
  std::vector<long long> flat_sorted(nkeys * KEY_SIZE, -1);
  if (rank == 0) {
    for (int i = 0; i < nkeys; i++) {
      for (int j = 0; j < KEY_SIZE; j++) flat_sorted[i * KEY_SIZE + j] = sorted_keys[i][j];
    }
  }
  MPI_Bcast(flat_sorted.data(), (int)flat_sorted.size(), MPI_LONG_LONG, 0, mesh.comm);
 
  // Reconstruct sorted_keys on all ranks
  if (rank != 0) {
    sorted_keys.resize(nkeys);
    for (int i = 0; i < nkeys; i++) {
      std::array<long long, KEY_SIZE> k;
      for (int j = 0; j < KEY_SIZE; j++) k[j] = flat_sorted[i * KEY_SIZE + j];
      sorted_keys[i] = k;
    }
  }
 
  // Assign deterministic element numbering (both REF and SUBMESH)
  SF::vector<mesh_int_t> & nbr_ref = mesh.register_numbering(SF::NBR_ELEM_REF);
  SF::vector<mesh_int_t> & nbr_sub = mesh.register_numbering(SF::NBR_ELEM_SUBMESH);
  nbr_ref.resize(mesh.l_numelem);
  nbr_sub.resize(mesh.l_numelem);
  for (size_t i = 0; i < mesh.l_numelem; i++) {
    auto k = make_key(i);
    auto it = std::lower_bound(sorted_keys.begin(), sorted_keys.end(), k);
    if (it == sorted_keys.end() || *it != k) {
      log_msg(0, 5, 0, "deterministic numbering failed to find key (rank %d, elem %zu)", rank, i);
      EXIT(1);
    }
    mesh_int_t gid = (mesh_int_t)(it - sorted_keys.begin());
    nbr_ref[i] = gid;
    nbr_sub[i] = gid;
  }
}
 
void EMI::setup_output()
{
  bool write_binary = false;
  std::string output_base = get_basename(param_globals::meshname);
  set_dir(OUTPUT);
 
  // write entire mesh
  sf_mesh & mesh = get_mesh(emi_msh);
  std::string output_file = output_base + "_e";
  log_msg(0, 0, 0, "Writing \"%s\" mesh: %s", mesh.name.c_str(), output_file.c_str());
  write_mesh_parallel(mesh, write_binary, output_file.c_str());
  // register output for overall phi on the entire mesh
  output_manager.register_output(parab_solver.ui, emi_msh, 1, param_globals::phiefile, "mV", NULL);
 
  // write membrane mesh
  sf_mesh & mesh_m = get_mesh(emi_surface_unique_face_msh);
  mesh_m.name = "Membrane";
  output_file = output_base + "_m";
  log_msg(0, 0, 0, "Writing \"%s\" mesh: %s", mesh_m.name.c_str(), output_file.c_str());
  write_mesh_parallel(mesh_m, write_binary, output_file.c_str());
  // register output for Vm on membrane interface
  // ensure deterministic element ordering in output
  if(get_permutation(emi_surface_unique_face_msh, ELEM_PETSC_TO_CANONICAL, 1) == NULL) {
    register_permutation(emi_surface_unique_face_msh, ELEM_PETSC_TO_CANONICAL, 1);
  }
  output_manager.register_output(parab_solver.vb_unique_face, emi_surface_unique_face_msh, 1, param_globals::vofile, "mV", NULL, true);
 
  if(param_globals::num_trace) {
    sf_mesh & imesh = get_mesh(emi_msh);
    open_trace(ion.miif, param_globals::num_trace, param_globals::trace_node, NULL, &imesh);
  }
 
  // initialize generic logger for IO timings per time_dt
  IO_stats.init_logger("IO_stats.dat");
}
 
void EMI::dump_matrices()
{
  std::string bsname = param_globals::dump_basename;
  std::string fn;
 
  set_dir(OUTPUT);
 
  fn = bsname + "_lhs.bin";
  parab_solver.lhs_emi->write(fn.c_str());
 
  fn = bsname + "_K.bin";
  parab_solver.stiffness_emi->write(fn.c_str());
 
  fn = bsname + "_B.bin";
  parab_solver.B->write(fn.c_str());
 
  fn = bsname + "_Bi.bin";
  parab_solver.Bi->write(fn.c_str());
 
  fn = bsname + "_BsM.bin";
  parab_solver.BsM->write(fn.c_str());
 
  fn = bsname + "_M.bin";
  parab_solver.mass_emi->write(fn.c_str());
 
  fn = bsname + "_Ms.bin";
  parab_solver.mass_surf_emi->write(fn.c_str());
 
}
 
double EMI::timer_val(const int timer_id)
{
  // determine
  int sidx = stimidx_from_timeridx(stimuli, timer_id);
  double val = 0.0;
  if(sidx != -1) {
    stimuli[sidx].value(val);
  }
  else
    val = std::nan("NaN");
 
  return val;
}
 
std::string EMI::timer_unit(const int timer_id)
{
  int sidx = stimidx_from_timeridx(stimuli, timer_id);
  std::string s_unit;
 
  if(sidx != -1)
    // found a timer-linked stimulus
    s_unit = stimuli[sidx].pulse.wave.f_unit;
 
  return s_unit;
}
 
void EMI::setup_solvers()
{
  set_dir(OUTPUT);
  parab_solver.init();
  parab_solver.rebuild_matrices(mtype_vol, *ion.miif, stimuli, logger);
 
  if(param_globals::dump2MatLab)
    dump_matrices();
}
 
void extract_unique_tag(SF::vector<mesh_int_t>& unique_tags)
{
  MPI_Comm comm = SF_COMM;
  int      size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  binary_sort(unique_tags);
  unique_resize(unique_tags); // unique_resize is done locally(for each rank)
  make_global(unique_tags, comm);
  binary_sort(unique_tags);
  unique_resize(unique_tags);
}
 
void compute_tags_per_rank(int num_tags, SF::vector<mesh_int_t>& num_tags_per_rank)
{
  MPI_Comm comm = SF_COMM;
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
  divide(num_tags, size, num_tags_per_rank);
}
 
void EMI::setup_EMI_mesh() 
{
  log_msg(0,0,0, "\n    *** Processing EMI mesh ***\n");
 
  const std::string basename = param_globals::meshname;
  const int verb = param_globals::output_level;
  std::map<mesh_t, sf_mesh> & mesh_registry = user_globals::mesh_reg;
  assert(mesh_registry.count(emi_msh) == 1); 
 
  set_dir(INPUT);
 
  sf_mesh & emi_mesh = mesh_registry[emi_msh];
  sf_mesh & emi_surfmesh_one_side = mesh_registry[emi_surface_msh];
  sf_mesh & emi_surfmesh_w_counter_face = mesh_registry[emi_surface_counter_msh];
  sf_mesh & emi_surfmesh_unique_face = mesh_registry[emi_surface_unique_face_msh];
 
  MPI_Comm comm = emi_mesh.comm;
 
  int size, rank;
  double t1, t2, s1, s2;
  const double total_setup_t0 = MPI_Wtime();
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
  //-----------------------------------------------------------------
  // Validate mesh dimension: EMI only supports 3D volumetric meshes
  //-----------------------------------------------------------------
  if (emi_mesh.l_numelem > 0) {
    // Check first local element type (all elements should be same type in CARP)
    SF::elem_t first_elem = emi_mesh.type[0];
    
    if (first_elem == SF::Line || first_elem == SF::Tri || first_elem == SF::Quad) {
      const char* type_name = (first_elem == SF::Line) ? "1D (Line)" :
                             (first_elem == SF::Tri)  ? "2D (Tri)" :
                                                         "2D (Quad)";
      if (rank == 0) {
        log_msg(0, 5, 0, "\n*** ERROR: EMI model requires a 3D volumetric mesh!");
        log_msg(0, 5, 0, "*** Current mesh element type: %s", type_name);
        log_msg(0, 5, 0, "*** EMI only supports 3D element types: Tetra, Pyramid, Prism, Hexa");
        log_msg(0, 5, 0, "*** Please provide a 3D mesh with volume elements.\n");
      }
      EXIT(EXIT_FAILURE);
    }
  }
 
  //-----------------------------------------------------------------
  // Step 1: READ *.intra and *.extra 
  //-----------------------------------------------------------------
  int total_num_tags = 0;
  if(verb) log_msg(NULL, 0, 0,"\nReading tags for extra and intra regions from input files");
  t1 = MPI_Wtime();
  {
    SF::vector<mesh_int_t> unique_extra_tags;
    SF::vector<mesh_int_t> unique_intra_tags;
 
    if(verb) log_msg(NULL, 0, 0,"Read extracellular tags");
    read_indices_global(unique_extra_tags,basename+".extra", comm);
    for(mesh_int_t tag:unique_extra_tags){
      parab_solver.extra_tags.insert(tag);
    }
 
    if(verb) log_msg(NULL, 0, 0,"Read intracellular tags");
    read_indices_global(unique_intra_tags,basename+".intra", comm);
    for(mesh_int_t tag:unique_intra_tags){
      parab_solver.intra_tags.insert(tag);
    }
 
    total_num_tags = parab_solver.extra_tags.size() + parab_solver.intra_tags.size();
    if(total_num_tags < size){
      log_msg(0,5,0, "\nThe number of unique tags on EMI mesh is smaller than number of processors!");
      EXIT(EXIT_FAILURE);
    }
    if(verb) log_msg(NULL, 0, 0,"\nextra_tags=%lu, intra_tags=%lu", parab_solver.extra_tags.size(), parab_solver.intra_tags.size());
  }
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 2: READ points from *.pts
  //-----------------------------------------------------------------
  SF::vector<mesh_real_t> pts;
  SF::vector<mesh_int_t>  ptsidx;
  SF::vector<mesh_int_t>  ptsData;
  if(verb) log_msg(NULL, 0, 0,"\nReading points with data on each vertex");
  t1 = MPI_Wtime();
  SF::read_points(basename, comm, pts, ptsidx);
  ptsData.resize(ptsidx.size());
  assert(ptsidx.size()==ptsData.size());
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  std::list< sf_mesh* > meshlist;
  meshlist.push_back(&emi_mesh);
  
  //-----------------------------------------------------------------
  // Step 3: Distrubute mesh based on tag or *.part
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0,"\nDistribute mesh based on tags");
  // should be replaced by scotch/pt-scotch for efficiency
  t1 = MPI_Wtime(); 
  distribute_elements_based_tags(emi_mesh, total_num_tags);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 4: insert coordinates to the mesh
  //-----------------------------------------------------------------
  // insert points
  if(verb) log_msg(NULL, 0, 0, "\nInserting points");
  t1 = MPI_Wtime();
  SF::insert_points(pts, ptsidx, meshlist);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 5: extract ptsdata based on the location of each vertex
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0, "\nCompute location of all dofs on emi mesh(inner dofs, memebarne, gap junction) saved into ptsData from original mesh");
  t1 = MPI_Wtime();
  compute_ptsdata_from_original_mesh( emi_mesh, 
                                      SF::NBR_REF,
                                      vertex2ptsdata, 
                                      parab_solver.extra_tags,
                                      parab_solver.intra_tags);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 6: extract faces and counterfaces on the interfaces (mem/gap) 
  // and mark the unique faces on surface mesh
  // generate required mapping between surface mesh and unique faces
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "\nExtract EMI surface mesh");
  hashmap::unordered_map<mesh_int_t, SF::emi_index_rank<mesh_int_t>> unused_map_elem_oneface_to_elem_uniqueFace;
  extract_face_based_tags(emi_mesh, SF::NBR_REF, vertex2ptsdata,
                          parab_solver.line_face, 
                          parab_solver.tri_face, 
                          parab_solver.quad_face, 
                          parab_solver.extra_tags, 
                          parab_solver.intra_tags, 
                          emi_surfmesh_one_side, emi_surfmesh_w_counter_face, emi_surfmesh_unique_face,
                          parab_solver.map_elem_uniqueFace_to_elem_oneface,
                          unused_map_elem_oneface_to_elem_uniqueFace);
  meshlist.push_back(&emi_surfmesh_one_side);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
  //-----------------------------------------------------------------
  // Step 7: register surface meshes based on the {typeof}_face with counter face
  //-----------------------------------------------------------------
  compute_surface_mesh_with_counter_face(emi_surfmesh_w_counter_face, SF::NBR_REF,
                                         parab_solver.line_face, 
                                         parab_solver.tri_face, 
                                         parab_solver.quad_face);
 
  compute_surface_mesh_with_unique_face(emi_surfmesh_unique_face, SF::NBR_REF,
                                         parab_solver.line_face, 
                                         parab_solver.tri_face, 
                                         parab_solver.quad_face,
                                         parab_solver.map_elem_uniqueFace_to_tags);
 
  //-----------------------------------------------------------------
  // Step 8: create a map between emi_surfmesh_w_counter_face and emi_surfmesh (both -> one)
  //-----------------------------------------------------------------
  SF::create_reverse_elem_mapping_between_surface_meshes(parab_solver.line_face,
                          parab_solver.tri_face,
                          parab_solver.quad_face,
                          parab_solver.vec_both_to_one_face,
                          comm);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 9: global number of interfaces only on emi_surfmesh without counter face
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0, "\ncompute global number of interface");
  size_t global_count_surf = 0;
  size_t numelem_surface = emi_surfmesh_one_side.l_numelem;
  unsigned long long local_count_surf = static_cast<unsigned long long>(numelem_surface);
  unsigned long long global_count_surf_ull = 0;
  MPI_Reduce(&local_count_surf, &global_count_surf_ull, 1, MPI_UNSIGNED_LONG_LONG, MPI_SUM, 0, MPI_COMM_WORLD);
  global_count_surf = static_cast<size_t>(global_count_surf_ull);
  if(verb && rank==0) fprintf(stdout, "global number of interfaces = %zu\n", global_count_surf);
 
  //-----------------------------------------------------------------
  // Step 10: submesh_numbering on the emi_mesh and generate parallel layout 
  //-----------------------------------------------------------------
  {
    SF::submesh_numbering<mesh_int_t,mesh_real_t> sub_numbering;
    sub_numbering(emi_mesh);
    emi_mesh.generate_par_layout();
  }
 
  SF::meshdata<mesh_int_t, mesh_real_t>     tmesh_backup_old = emi_mesh;
 
  //-----------------------------------------------------------------
  // Step 11: make a map bwteen key(vertex,tag)-> value(dof) after decoupling interfaces defined on emi mesh
  //-----------------------------------------------------------------
  // During interface decoupling and the introduction of new degrees of freedom (DoFs),
  // we assume that the mesh partitioning is performed at least on a per-tag basis.
  // This means that all elements sharing the same tag number belong to the same rank.
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "\ndecouple emi interfaces");
  if(verb) log_msg(NULL, 0, 0, "\tcompute map oldIdx to dof");
  compute_map_vertex_to_dof(emi_mesh, SF::NBR_REF, vertex2ptsdata, parab_solver.extra_tags, parab_solver.map_vertex_tag_to_dof);
  
  //-----------------------------------------------------------------
  // Step 12: complete the map bwteen key(vertex,tag)-> value(dof) for counter faces defined on emi mesh
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0, "\tcomplete map oldIdx to dof with counter interface");
  // add to the map the counter part of the interface
  complete_map_vertex_to_dof_with_counter_face(parab_solver.line_face, parab_solver.tri_face, parab_solver.quad_face, parab_solver.map_vertex_tag_to_dof);
  
  //-----------------------------------------------------------------
  // Step 13: update emi mesh with new dofs, so decoupling the interfaces applied on emi mesh
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0, "\tupdate mesh with dof");
  update_emi_mesh_with_dofs(emi_mesh, SF::NBR_REF, parab_solver.map_vertex_tag_to_dof, parab_solver.dof2vertex);  
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 14: set map_vertex_tag_to_dof_petsc where the key(vertex,tag) and value(dof,petsc) 
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "\nInitialize petsc =0 for  map<oldIdx,tag> -><dof, petsc>");
  // Iterate over map to assign the (oldIndx, tag) -> (dof, petsc) atm petsc = 0
  for(const auto & key_value : parab_solver.map_vertex_tag_to_dof) 
  {
    mesh_int_t gIndex_old  = key_value.first.first;
    mesh_int_t tag_old  = key_value.first.second;
    mesh_int_t dof = key_value.second;
 
    std::pair <mesh_int_t,mesh_int_t> dof_petsc = std::make_pair(dof,-1); // petsc numbering ...
    parab_solver.map_vertex_tag_to_dof_petsc.insert({key_value.first,dof_petsc});
  }
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 16: insert the coordinates on new dofs to the emi mesh
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Inserting points and ptsData of dofs to emi_mesh");
  insert_points_ptsData_to_dof(tmesh_backup_old, emi_mesh, SF::NBR_REF, parab_solver.dof2vertex, vertex2ptsdata, dof2ptsData);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 17: submesh_numbering on the emi_mesh and generate parallel layout after decoupling the interfaces on the emi mesh 
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Generating unique PETSc numberings");
  {
    SF::submesh_numbering<mesh_int_t,mesh_real_t> sub_numbering;
    sub_numbering(emi_mesh);
    emi_mesh.generate_par_layout(); 
  }
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
  
  //-----------------------------------------------------------------
  // Step 18: register petsc_numbering on emi mesh
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Generating unique PETSc numberings");
  {
    SF::petsc_numbering<mesh_int_t,mesh_real_t> petsc_numbering(emi_mesh.pl);
    petsc_numbering(emi_mesh);
  }
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 19: complete map_vertex_tag_to_dof_petsc where the key(vertex,tag) and value(dof,petsc) for counter faces on emi mesh
  //-----------------------------------------------------------------
  if(verb) log_msg(NULL, 0, 0, "Updating the map between indices to PETSc numberings");
  t1 = MPI_Wtime();
  update_map_indices_to_petsc(emi_mesh, SF::NBR_REF, SF::NBR_PETSC, parab_solver.extra_tags, parab_solver.map_vertex_tag_to_dof_petsc, parab_solver.dof2vertex, parab_solver.elemTag_emi_mesh);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 20: update value of map {type_of}_face with new dofs after decoupling
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Updating surface mesh with dof");
  update_faces_on_surface_mesh_after_decoupling_with_dofs(emi_surfmesh_one_side, parab_solver.map_vertex_tag_to_dof, parab_solver.line_face, parab_solver.tri_face, parab_solver.quad_face);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 21: register emi_surfmesh for SF::NBR_ELEM_REF, NBR_REF, SF::NBR_SUBMESH
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Layout for element of EMI surfmesh");
  SF::vector<mesh_int_t> & emi_surfmesh_elem =  emi_surfmesh_one_side.register_numbering(SF::NBR_ELEM_REF);
  SF::vector<long int> layout;
  SF::layout_from_count<long int>(emi_surfmesh_one_side.l_numelem, layout, emi_surfmesh_one_side.comm);
  size_t count = layout[rank+1] - layout[rank];
  emi_surfmesh_elem.resize(count);
  for (int i = 0; i < count; ++i){
    emi_surfmesh_elem[i] = layout[rank]+i;
  }
 
  emi_surfmesh_one_side.localize(SF::NBR_REF);
  emi_surfmesh_one_side.register_numbering(SF::NBR_SUBMESH);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 22: register emi_surfmesh_w_counter_face for SF::NBR_ELEM_REF, NBR_REF, SF::NBR_SUBMESH
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Layout for element of EMI surfmesh w counter face");
  SF::vector<mesh_int_t> & emi_surfmesh_counter_elem =  emi_surfmesh_w_counter_face.register_numbering(SF::NBR_ELEM_REF);
  SF::vector<long int> layout_counter;
  SF::layout_from_count<long int>(emi_surfmesh_w_counter_face.l_numelem, layout_counter, emi_surfmesh_w_counter_face.comm);
  size_t count_counter = layout_counter[rank+1] - layout_counter[rank];
  emi_surfmesh_counter_elem.resize(count_counter);
  for (int i = 0; i < count_counter; ++i){
    emi_surfmesh_counter_elem[i] = layout_counter[rank]+i;
  }
  emi_surfmesh_w_counter_face.localize(SF::NBR_REF);
  emi_surfmesh_w_counter_face.register_numbering(SF::NBR_SUBMESH);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
 
  //-----------------------------------------------------------------
  // Step 23: register emi_surfmesh_w_counter_face for SF::NBR_ELEM_REF, NBR_REF, SF::NBR_SUBMESH
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Layout for element of EMI surfmesh w counter face");
  SF::vector<mesh_int_t> & emi_surfmesh_unique_elem =  emi_surfmesh_unique_face.register_numbering(SF::NBR_ELEM_REF);
  SF::vector<long int> layout_unique;
  SF::layout_from_count<long int>(emi_surfmesh_unique_face.l_numelem, layout_unique, emi_surfmesh_unique_face.comm);
  size_t count_unique = layout_unique[rank+1] - layout_unique[rank];
  emi_surfmesh_unique_elem.resize(count_unique);
  for (int i = 0; i < count_unique; ++i){
    emi_surfmesh_unique_elem[i] = layout_unique[rank]+i;
  }
  emi_surfmesh_unique_face.localize(SF::NBR_REF);
  emi_surfmesh_unique_face.register_numbering(SF::NBR_SUBMESH);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 24: insert coorinates to emi_surfmesh_one_side and emi_surfmesh_w_counter_face 
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Inserting points to EMI surfmesh");
  insert_points_to_surface_mesh(tmesh_backup_old, emi_surfmesh_one_side, SF::NBR_REF, parab_solver.dof2vertex, parab_solver.extra_tags, parab_solver.elemTag_surface_mesh);
  insert_points_to_surface_mesh(tmesh_backup_old, emi_surfmesh_w_counter_face, SF::NBR_REF, parab_solver.dof2vertex, parab_solver.extra_tags, parab_solver.elemTag_surface_w_counter_mesh); 
  insert_points_to_surface_mesh(tmesh_backup_old, emi_surfmesh_unique_face, SF::NBR_REF, parab_solver.dof2vertex);
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 25: submesh_numbering & petc numebring for both emi_surfmesh_one_side and emi_surfmesh_w_counter_face 
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Generating submesh_numbering and PETSc numberings for surface mesh");
  {
    SF::submesh_numbering<mesh_int_t,mesh_real_t> sub_numbering;
    sub_numbering(emi_surfmesh_one_side);
    emi_surfmesh_one_side.generate_par_layout(); 
 
    SF::petsc_numbering<mesh_int_t,mesh_real_t> petsc_numbering(emi_surfmesh_one_side.pl);
    petsc_numbering(emi_surfmesh_one_side);
  }
 
  {
    SF::submesh_numbering<mesh_int_t,mesh_real_t> sub_numbering;
    sub_numbering(emi_surfmesh_w_counter_face);
    emi_surfmesh_w_counter_face.generate_par_layout();
 
    SF::petsc_numbering<mesh_int_t,mesh_real_t> petsc_numbering(emi_surfmesh_w_counter_face.pl);
    petsc_numbering(emi_surfmesh_w_counter_face);
  }
 
  {
    SF::submesh_numbering<mesh_int_t,mesh_real_t> sub_numbering;
    sub_numbering(emi_surfmesh_unique_face);
    emi_surfmesh_unique_face.generate_par_layout();
 
    SF::petsc_numbering<mesh_int_t,mesh_real_t> petsc_numbering(emi_surfmesh_unique_face.pl);
    petsc_numbering(emi_surfmesh_unique_face);
    // Ensure deterministic element ordering across MPI for output
    assign_deterministic_elem_numbering(emi_surfmesh_unique_face);
  }
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  //-----------------------------------------------------------------
  // Step 26: assign petsc numbering to map_vertex_tag_to_dof_petsc where the key(vertex,tag) and value(dof,petsc) 
  //-----------------------------------------------------------------
  t1 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "assign PETSc numbering for new faces");
  SF::assign_petsc_on_counter_face(parab_solver.map_vertex_tag_to_dof_petsc,comm);
  {
    hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, std::pair<mesh_int_t,mesh_int_t>>::iterator it;
    for (it = parab_solver.map_vertex_tag_to_dof_petsc.begin(); it != parab_solver.map_vertex_tag_to_dof_petsc.end(); it++)
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old = it->first;
      std::pair <mesh_int_t,mesh_int_t> dof_petsc = it->second;
      parab_solver.dof2petsc[dof_petsc.first] = dof_petsc.second;
      parab_solver.petsc2dof[dof_petsc.second] = dof_petsc.first;
    }
  }
  
  // DEBUG: Check for invalid PETSc indices after exchange
  #ifdef EMI_DEBUG_MESH
  {
    int invalid_count = 0;
    for (const auto& [key, val] : parab_solver.map_vertex_tag_to_dof_petsc) {
      if (val.second < 0) {
        invalid_count++;
        if (invalid_count <= 3) {
          fprintf(stderr, "RANK %d INVALID: vertex=%ld tag=%ld dof=%ld petsc=%ld\n", 
                  rank, (long)key.first, (long)key.second, (long)val.first, (long)val.second);
        }
      }
    }
    fprintf(stderr, "RANK %d: After Step 26: %d invalid PETSc indices out of %zu total\n", 
            rank, invalid_count, parab_solver.map_vertex_tag_to_dof_petsc.size());
    fflush(stderr);
  }
  #endif
  
  t2 = MPI_Wtime();
  if(verb) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
  
  //-----------------------------------------------------------------
  // Step 27: update counter faces to the map{type_of}_face
  //-----------------------------------------------------------------
  added_counter_faces_to_map(parab_solver.line_face, parab_solver.tri_face, parab_solver.quad_face);
  const double total_setup = MPI_Wtime() - total_setup_t0;
  log_msg(0,0,0, "Total setup_EMI_mesh processing time: %.5f sec.", float(total_setup));
  
  log_msg(0,0,0, "\n    *** EMI mesh processing Done ***\n");
}
 
void distribute_elements_based_tags(SF::meshdata<mesh_int_t, mesh_real_t>& mesh,
                                    int total_num_tags)
{
  MPI_Comm comm = mesh.comm;
  int size, rank;
  double t1, t2, s1, s2;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
  const int verb = param_globals::output_level;
 
  if(total_num_tags < size)
  {
    PetscPrintf(PETSC_COMM_WORLD,"\nThe number of processors should be less than the number of tags, size = %d & ntags = %d !!!\n",
                size, total_num_tags);
    cleanup_and_exit();
  }
 
  if(verb==10) log_msg(NULL, 0, 0,"\ncompute the number of tags which belongs to one rank");
  // compute the number of tags per rank
  SF::vector<mesh_int_t> ntags_per_rank;
  t1 = MPI_Wtime();
  compute_tags_per_rank(total_num_tags, ntags_per_rank); 
  t2 = MPI_Wtime();
  if(verb==10) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
  // generate global destination w.r.t. tags
  // redistribute elements based on new partition 
  t1 = MPI_Wtime(); 
  SF::vector<mesh_int_t> part_based_Tags(mesh.l_numelem);
  partition_based_tags(total_num_tags, mesh.tag, ntags_per_rank, part_based_Tags);
  SF::redistribute_elements(mesh,part_based_Tags);
  // permute elements locally first based on the tag then element index
  permute_mesh_locally_based_on_tag_elemIdx(mesh);
  t2 = MPI_Wtime(); 
  if(verb==10) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
}
 
// In the EMI model we cannot have partition boundaries crossing cells.
// Therefore this function assigns a partition to each tag. All elements with this tag will go to the same partition.
void partition_based_tags(int num_tags,
                          SF::vector<mesh_int_t> tag,
                          SF::vector<mesh_int_t> num_tags_per_rank,
                          SF::vector<mesh_int_t>& part)
{
  const int verb = param_globals::output_level;
  MPI_Comm comm = SF_COMM;
  int      size, rank;
  MPI_Comm_size(comm, &size);
  MPI_Comm_rank(comm, &rank);
 
  // we need to have tags_to_rank as a map globally.
  hashmap::unordered_map<mesh_int_t, mesh_int_t> tags_to_rank_map; 
  // Try to load the mapping from a .part file and fall back on the internal mapping if there is no such file.
  // If the mapping is read from file, the num_tags_per_rank argument is not used.
  if (!load_partitions_from_file(tags_to_rank_map, num_tags, comm)) {
    if(verb==10) log_msg(NULL, 0, 0,"\ncompute the number of unique tags");
    SF::vector<mesh_int_t> unique_tags = tag;
    double t1 = MPI_Wtime();
    extract_unique_tag(unique_tags);
    double t2 = MPI_Wtime();
    if(verb==10) log_msg(NULL, 0, 0, "Done in %f sec.", float(t2 - t1));
 
    if (unique_tags.size() != static_cast<size_t>(num_tags)) {
      log_msg(0,5,0, "\nerror: the number of tags in the EMI mesh does not match the total tags from *.extra and *.intra");
      EXIT(1);
    }
    map_tags_to_rank(size, unique_tags, num_tags_per_rank, tags_to_rank_map);
  }
 
  for (size_t i = 0; i < part.size(); ++i) {
    if(tags_to_rank_map.count(tag[i]))
      part[i] = tags_to_rank_map[tag[i]];
  }
}
 
void map_tags_to_rank(int size, const SF::vector<mesh_int_t> & unique_tags, const SF::vector<mesh_int_t> & num_tags_per_rank, hashmap::unordered_map<mesh_int_t, mesh_int_t> &tags_to_rank_map)
{
  // the goal is to map the tag as a key -> rank, should be global one
  mesh_int_t t = 0;  // tag index
  for (size_t r = 0; r < size; ++r) {
    for (size_t count = 0; count < num_tags_per_rank[r];) {
      mesh_int_t tag =  unique_tags[t];
 
      // // here we consider each tag belongs to only one rank.
      tags_to_rank_map.insert({tag, r});  
      count +=1; 
      t+=1;
    }
  }
}
 
// This function creates a tag-to-rank mapping based only on the numerical order of the tags.
bool load_partitions_from_file(hashmap::unordered_map<mesh_int_t, mesh_int_t>& tags_to_rank_map,
                               int expected_num_tags,
                               MPI_Comm comm)
{
  FILE_SPEC tags_to_part = nullptr;
  
  SF::vector<int> parts;
  const std::string basename = param_globals::meshname;
  FILE* fd = fopen((basename + ".part").c_str(), "r");
  if (fd != NULL) { // check if the file exists otherwise read_indices_global() crash
    read_indices_global(parts, basename + ".part", comm);
    fclose(fd);
  } 
  
  if (parts.size() == 0) return false;  // file does not exist or is empty
 
  if (parts.size() % 2 != 0) {
    log_msg(0,5,0, "\nThe part file should contain 2 lines per tag, one for the tag number and the next for its associated partition number.!");
    EXIT(1);
  }
 
  for(int i = 0; i < parts.size(); i+=2) {
    tags_to_rank_map.insert({parts[i], parts[i + 1]});
  }
 
  if (tags_to_rank_map.size() != static_cast<size_t>(expected_num_tags)) {
    log_msg(0,5,0, "\nerror: the number of tags in the .part file does not match the total tags from *.extra and *.intra");
    EXIT(1);
  }
 
  return true;
}
 
void permute_mesh_locally_based_on_tag_elemIdx(SF::meshdata<mesh_int_t, mesh_real_t>& mesh)
{
  mesh.globalize(SF::NBR_REF);
  SF::meshdata<mesh_int_t, mesh_real_t>     tmesh = mesh;
  SF::vector<mesh_int_t> perm;
  interval(perm, 0, mesh.tag.size());
 
  SF::vector<mesh_int_t>  tags    = tmesh.tag;
  SF::vector<mesh_int_t>& elemIdx = tmesh.get_numbering(SF::NBR_ELEM_REF);
  binary_sort_sort_copy(tags, elemIdx, perm);
  permute_mesh(tmesh, mesh, perm);
 
  mesh.localize(SF::NBR_REF);
}
 
}  // namespace opencarp
 
#endif