electrics_eikonal.cc

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
// ----------------------------------------------------------------------------

#include "electrics_eikonal.h"
#include "basics.h"
#include "petsc_utils.h"
#include "timers.h"
#include "stimulate.h"

#include "SF_init.h" 

namespace opencarp
{

void Eikonal::initialize()
{
  double t1, t2;
  get_time(t1);

  set_dir(OUTPUT);

  // open logger
  logger = f_open("eikonal.log", param_globals::experiment != 4 ? "w" : "r");

  // setup mappings between extra and intra grids, algebraic and nodal,
  // and between PETSc and canonical orderings
  setup_mappings();

  eik_tech = static_cast<Eikonal::eikonal_t>(param_globals::dream.solve);

  // the ionic physics is currently triggered from inside the Electrics to have tighter
  // control over it.
  switch (eik_tech) {
    case EIKONAL: break;
    default:
      ion.logger = logger;
      ion.initialize();
  }

  // set up Intracellular tissue
  set_elec_tissue_properties(mtype, intra_grid, logger);
  region_mask(intra_elec_msh, mtype[intra_grid].regions, mtype[intra_grid].regionIDs, true, "gregion_i");

  // 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");

  // electrics stimuli setup
  setup_stimuli();

  // set up the linear equation systems. this needs to happen after the stimuli have been
  // set up, since we need boundary condition info
  setup_solvers();

  // the next setup steps require the solvers to be set up, since they use the matrices
  // generated by those

  // initialize the LATs detector
  switch (eik_tech) {
    case EIKONAL: break; // not available for pure eikonal solve
    default:
      lat.init(*parab_solver.Vmv, *phie_dummy, 0, eikonal_phys);
  }

  // prepare the electrics output. we skip it if we do post-processing
  if (param_globals::experiment != EXP_POSTPROCESS)
    setup_output();

  if (param_globals::prepacing_bcl > 0)
    prepace();

  this->initialize_time += timing(t2, t1);
}

void Eikonal::set_elec_tissue_properties(MaterialType* mtype, Eikonal::grid_t g, FILE_SPEC logger)
{
  MaterialType* m = mtype + g;

  // initialize random conductivity fluctuation structure with PrM values
  m->regions.resize(param_globals::num_gregions);

  const char* grid_name = g == Eikonal::intra_grid ? "intracellular" : "extracellular";
  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();

  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);
    }

    reg->regname  = strdup(param_globals::gregion[i].name);
    reg->regID    = i;
    reg->nsubregs = param_globals::gregion[i].num_IDs;
    if (!reg->nsubregs)
      reg->subregtags = NULL;
    else {
      reg->subregtags = new int[reg->nsubregs];

      for (int j = 0; j < reg->nsubregs; j++)
        reg->subregtags[j] = param_globals::gregion[i].ID[j];
    }

    // describe material in given region
    elecMaterial* emat  = new elecMaterial();
    emat->material_type = ElecMat;

    emat->InVal[0] = param_globals::gregion[i].g_il;
    emat->InVal[1] = param_globals::gregion[i].g_it;
    emat->InVal[2] = param_globals::gregion[i].g_in;

    emat->ExVal[0] = param_globals::gregion[i].g_el;
    emat->ExVal[1] = param_globals::gregion[i].g_et;
    emat->ExVal[2] = param_globals::gregion[i].g_en;

    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 ((g == Eikonal::intra_grid && strlen(param_globals::gi_scale_vec)) ||
      (g == Eikonal::extra_grid && strlen(param_globals::ge_scale_vec))) {
    mesh_t      mt   = g == Eikonal::intra_grid ? intra_elec_msh : extra_elec_msh;
    sf_mesh&    mesh = get_mesh(mt);
    const char* file = g == Eikonal::intra_grid ? param_globals::gi_scale_vec : param_globals::ge_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(g == Eikonal::intra_grid ? param_globals::gi_scale_vec : param_globals::ge_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 Eikonal::setup_mappings()
{
  bool intra_exits = mesh_is_registered(intra_elec_msh), extra_exists = mesh_is_registered(extra_elec_msh);
  assert(intra_exits);
  const int dpn = 1;

  // It may be that another physic (e.g. ionic models) has already computed the intracellular mappings,
  // thus we first test their existence
  if (get_scattering(intra_elec_msh, ALG_TO_NODAL, dpn) == NULL) {
    log_msg(logger, 0, 0, "%s: Setting up intracellular algebraic-to-nodal scattering.", __func__);
    register_scattering(intra_elec_msh, ALG_TO_NODAL, dpn);
  }
  if (get_permutation(intra_elec_msh, PETSC_TO_CANONICAL, dpn) == NULL) {
    log_msg(logger, 0, 0, "%s: Setting up intracellular PETSc to canonical permutation.", __func__);
    register_permutation(intra_elec_msh, PETSC_TO_CANONICAL, dpn);
  }
}

void Eikonal::compute_step()
{
  double t1, t2;
  get_time(t1);

  switch (eik_tech) {
    case EIKONAL: solve_EIKONAL(); break;
    case DREAM: solve_DREAM(); break;
    default: solve_RE(); break;
  }

  if (user_globals::tm_manager->trigger(iotm_console) && eik_tech != EIKONAL) {
    // 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) && eik_tech != EIKONAL)
    dump_trace(ion.miif, user_globals::tm_manager->time);
}

void Eikonal::output_step()
{
  double t1, t2;
  get_time(t1);

  switch (eik_tech) {
    case EIKONAL: break;
    default:
      assemble_sv_gvec(gvec, ion.miif);
      output_manager_time.write_data();  // does not exist in EIKONAL
  }

  if (do_output_eikonal) {
    output_manager_cycle.write_data();
    do_output_eikonal = false;
  }

  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 Eikonal::destroy()
{
  // close logger
  f_close(logger);

  switch (eik_tech) {
    case EIKONAL:
      output_manager_cycle.close_files_and_cleanup();
      break;
    default:
      // output LAT data
      lat.output_initial_activations();
      output_manager_time.close_files_and_cleanup();
      output_manager_cycle.close_files_and_cleanup();
      ion.destroy();
  }
}

void Eikonal::setup_stimuli()
{
  // initialize basic stim info data (used units, supported types, etc)
  init_stim_info();
  bool dumpTrace = true;

  if (dumpTrace) set_dir(OUTPUT);

  stimuli.resize(param_globals::num_stim);
  for (int i = 0; i < param_globals::num_stim; i++) {
    // construct new stimulus
    stimulus& s = stimuli[i];

    if (param_globals::stim[i].crct.type != 0 && param_globals::stim[i].crct.type != 9) {
      // In the Eikonal class, only the intracellular domain is registered.
      // Therefore, only I_tm and Vm_clmp are compatible.
      log_msg(NULL, 5, 0, "%s error: stimulus of type %i is incompatible with the eikonal model! Use I_tm or Vm_clmp instead. Aborting!", __func__, s.phys.type);
      EXIT(EXIT_FAILURE);
    }

    if (user_globals::using_legacy_stimuli)
      s.translate(i);

    s.setup(i);

    log_msg(NULL, 2, 0, "Only geometry, start time, npls, and bcl of stim[%i] are used", i);

    if (dumpTrace && get_rank() == 0)
      s.pulse.wave.write_trace(s.name + ".trc");
  }

  eik_solver.set_stimuli(stimuli);
}

void Eikonal::stimulate_intracellular()
{
  parabolic_solver& ps = parab_solver;

  // iterate over stimuli
  for (stimulus& s : stimuli) {
    if (s.is_active()) {
      // for active stimuli, deal with the stimuli-type specific stimulus application
      switch (s.phys.type) {
        case I_tm: {
          if (param_globals::operator_splitting) {
            apply_stim_to_vector(s, *ps.Vmv, true);
          } else {
            SF_real        Cm = 1.0;
            timer_manager& tm = *user_globals::tm_manager;
            SF_real        sc = tm.time_step / Cm;

            ps.Irhs->set(0.0);
            apply_stim_to_vector(s, *ps.Irhs, true);

            *ps.tmp_i1 = *ps.IIon;
            *ps.tmp_i1 -= *ps.Irhs;
            *ps.tmp_i1 *= sc;  // tmp_i1 = sc * (IIon - Irhs)

            // add ionic, transmembrane and intracellular currents to rhs
            if (param_globals::parab_solve != parabolic_solver::EXPLICIT)
              ps.mass_i->mult(*ps.tmp_i1, *ps.Irhs);
            else
              *ps.Irhs = *ps.tmp_i1;
          }
          break;
        }

        case Illum: {
          sf_vec* illum_vec = ion.miif->gdata[limpet::illum];

          if (illum_vec == NULL) {
            log_msg(0, 5, 0, "Cannot apply illumination stim: global vector not present!");
            EXIT(EXIT_FAILURE);
          } else {
            apply_stim_to_vector(s, *illum_vec, false);
          }

          break;
        }

        default: break;
      }
    }
  }
}

void Eikonal::clamp_Vm()
{
  for (stimulus& s : stimuli) {
    if (s.phys.type == Vm_clmp && s.is_active())
      apply_stim_to_vector(s, *parab_solver.Vmv, false);
  }
}

void Eikonal::setup_output()
{
  int                     rank    = get_rank();
  SF::vector<mesh_int_t>* restr_i = NULL;
  SF::vector<mesh_int_t>* restr_e = NULL;
  set_dir(OUTPUT);

  setup_dataout(param_globals::dataout_i, param_globals::dataout_i_vtx, intra_elec_msh,
                restr_i, param_globals::num_io_nodes > 0);

  if (param_globals::dataout_i) {
    switch (eik_tech) {
      case EIKONAL:
        output_manager_cycle.register_output(eik_solver.AT, intra_elec_msh, 1, param_globals::dream.output.atfile, "ms", restr_i);
        break;
      case DREAM:
        output_manager_time.register_output(parab_solver.Vmv, intra_elec_msh, 1, param_globals::vofile, "mV", restr_i);
        output_manager_cycle.register_output(eik_solver.AT, intra_elec_msh, 1, param_globals::dream.output.atfile, "ms", restr_i);
        output_manager_cycle.register_output(eik_solver.RT, intra_elec_msh, 1, param_globals::dream.output.rtfile, "ms", restr_i);
        if (strcmp(param_globals::dream.output.idifffile, "") != 0) {
          output_manager_time.register_output(eik_solver.Idiff, intra_elec_msh, 1, param_globals::dream.output.idifffile, "muA/cm2", restr_i);
        }
        break;
      default:
        output_manager_time.register_output(parab_solver.Vmv, intra_elec_msh, 1, param_globals::vofile, "mV", restr_i);
        output_manager_cycle.register_output(eik_solver.AT, intra_elec_msh, 1, param_globals::dream.output.atfile, "ms", restr_i);
        if (strcmp(param_globals::dream.output.idifffile, "") != 0) {
          output_manager_time.register_output(eik_solver.Idiff, intra_elec_msh, 1, param_globals::dream.output.idifffile, "muA/cm2", restr_i);
        }
    }
  }

  init_sv_gvec(gvec, ion.miif, *parab_solver.Vmv, output_manager_time);

  if (param_globals::num_trace) {
    sf_mesh& imesh = get_mesh(intra_elec_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 Eikonal::dump_matrices()
{
  std::string bsname = param_globals::dump_basename;
  std::string fn;

  set_dir(OUTPUT);

  // dump monodomain matrices
  if (param_globals::parab_solve == 1) {
    // using Crank-Nicolson
    fn = bsname + "_Ki_CN.bin";
    parab_solver.lhs_parab->write(fn.c_str());
  }
  fn = bsname + "_Ki.bin";
  parab_solver.rhs_parab->write(fn.c_str());

  fn = bsname + "_Mi.bin";
  parab_solver.mass_i->write(fn.c_str());
}

double Eikonal::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 Eikonal::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 Eikonal::setup_solvers()
{
  set_dir(OUTPUT);

  switch (eik_tech) {
    case EIKONAL:
      eik_solver.init();
      eik_solver.set_initial_cv();
      break;
    default:
      parab_solver.init();
      parab_solver.rebuild_matrices(mtype, *ion.miif, logger);
      eik_solver.init();
      eik_solver.set_initial_cv();
      if (param_globals::dump2MatLab) dump_matrices();
  }
}

void Eikonal::checkpointing()
{
  const timer_manager& tm = *user_globals::tm_manager;

  // regular user selected state save
  if (tm.trigger(iotm_chkpt_list)) {
    char        save_fnm[1024];
    const char* tsav_ext = get_tsav_ext(tm.time);

    snprintf(save_fnm, sizeof save_fnm, "%s.%s.roe", param_globals::write_statef, tsav_ext);

    ion.miif->dump_state(save_fnm, tm.time, intra_elec_msh, false, GIT_COMMIT_COUNT);
    eik_solver.save_eikonal_state(tsav_ext);
  }

  // checkpointing based on interval
  if (tm.trigger(iotm_chkpt_intv)) {
    char save_fnm[1024];
    snprintf(save_fnm, sizeof save_fnm, "checkpoint.%.1f.roe", tm.time);
    ion.miif->dump_state(save_fnm, tm.time, intra_elec_msh, false, GIT_COMMIT_COUNT);
  }
}

void Eikonal::prepace()
{
  // Default stimulation parameters for the single cells
  // To be eventually be replaced by user input or info from bench
  float stimdur = 1.;
  float stimstr = 60.;

  log_msg(NULL, 0, 0, "Using activation times from file %s to distribute prepacing states\n",
          param_globals::prepacing_lats);
  log_msg(NULL, 1, 0, "Assuming stimulus strength %f uA/uF with duration %f ms for prepacing\n",
          stimstr, stimdur);

  Eikonal*          elec = (Eikonal*)get_physics(elec_phys);
  limpet::MULTI_IF* miif = elec->ion.miif;

  const sf_mesh& mesh = get_mesh(intra_elec_msh);
  sf_vec*        read_lats;
  SF::init_vector(&read_lats, mesh, 1, sf_vec::algebraic);

  // read in the global distributed vector of all activation times
  set_dir(INPUT);
  size_t numread = read_lats->read_ascii(param_globals::prepacing_lats);
  if (numread == 0) {
    log_msg(NULL, 5, 0, "Failed reading required LATs! Skipping prepacing!");
    return;
  }
  set_dir(OUTPUT);

  SF::scattering* sc = get_permutation(intra_elec_msh, PETSC_TO_CANONICAL, 1);
  assert(sc != NULL);

  // permute in-place to petsc permutation
  bool forward = false;
  (*sc)(*read_lats, forward);

  // take care of negative LAT values
  {
    SF_real* lp = read_lats->ptr();
    for (int i = 0; i < read_lats->lsize(); i++)
      if (lp[i] < 0.0) lp[i] = param_globals::tend + 10.0;

    read_lats->release_ptr(lp);
  }

  // make LATs relative and figure out the first LAT
  // so we know when to save state of each point
  SF_real LATmin = read_lats->min();

  if (LATmin < 0.0) {
    log_msg(0, 3, 0, "LAT data is not complete. Skipping prepacing.");
    return;
  }

  SF_real offset  = floor(LATmin / param_globals::prepacing_bcl) * param_globals::prepacing_bcl;
  SF_real last_tm = param_globals::prepacing_bcl * param_globals::prepacing_beats;

  // compute read_lats[i] = last_tm - (read_lats[i] - offset)
  *read_lats += -offset;
  *read_lats *= -1.;
  *read_lats += last_tm;

  miif->getRealData();
  SF_real* save_tm = read_lats->ptr();
  SF_real* vm      = miif->gdata[limpet::Vm]->ptr();

  for (int ii = 0; ii < miif->N_IIF; ii++) {
    if (!miif->N_Nodes[ii]) continue;

    // create sorted array of save times.
    SF::vector<SF::mixed_tuple<double, int>> sorted_save(miif->N_Nodes[ii]);  // v1 = time, v2 = index
    for (int kk = 0; kk < miif->N_Nodes[ii]; kk++) {
      sorted_save[kk].v1 = save_tm[miif->NodeLists[ii][kk]];
      sorted_save[kk].v2 = kk;
    }
    std::sort(sorted_save.begin(), sorted_save.end());

    size_t lastidx = sorted_save.size() - 1;
    int    paced   = sorted_save[lastidx].v2;  // IMP index of latest node
    int    csav    = 0;

    for (double t = 0; t < sorted_save[lastidx].v1; t += param_globals::dt) {
      if (fmod(t, param_globals::prepacing_bcl) < stimdur &&
          t < param_globals::prepacing_bcl * param_globals::prepacing_beats - 1)
        miif->ldata[ii][limpet::Vm][paced] += stimstr * param_globals::dt;

      compute_IIF(*miif->IIF[ii], miif->ldata[ii], paced);

      // Vm update always happens now outside of the imp
      miif->ldata[ii][limpet::Vm][paced] -= miif->ldata[ii][limpet::Iion][paced] * param_globals::dt;
      vm[miif->NodeLists[ii][paced]] = miif->ldata[ii][limpet::Vm][paced];

      while (csav < miif->N_Nodes[ii] - 1 && t >= sorted_save[csav].v1)
        dup_IMP_node_state(*miif->IIF[ii], paced, sorted_save[csav++].v2, miif->ldata[ii]);
    }

    // get nodes which may be tied for last
    while (csav < miif->N_Nodes[ii] - 1)
      dup_IMP_node_state(*miif->IIF[ii], paced, sorted_save[csav++].v2, miif->ldata[ii]);
    // ipdate global Vm vector
    for (int k = 0; k < miif->N_Nodes[ii]; k++) vm[miif->NodeLists[ii][k]] = miif->ldata[ii][limpet::Vm][k];
  }

  read_lats->release_ptr(save_tm);
  miif->gdata[limpet::Vm]->release_ptr(vm);
  miif->releaseRealData();
}

void Eikonal::solve_EIKONAL()
{
  if (user_globals::tm_manager->time == 0) {
    eik_solver.FIM();
    eik_solver.stats.log_stats(user_globals::tm_manager->time, false);
    eik_solver.nodeData.log_stats(user_globals::tm_manager->time, false);
    do_output_eikonal = true;
  }
}

void Eikonal::solve_RE()
{
  if (user_globals::tm_manager->time == 0) {
    eik_solver.FIM();
    eik_solver.stats.log_stats(user_globals::tm_manager->time, false);
    eik_solver.nodeData.log_stats(user_globals::tm_manager->time, false);
    do_output_eikonal = true;
  }
  solve_RD();
}

void Eikonal::solve_DREAM()
{
  if (eik_solver.determine_model_to_run(user_globals::tm_manager->time)) {
    solve_RD();
  } else {
    if (user_globals::tm_manager->time != 0) {
      eik_solver.update_repolarization_times(this->ion);
    }

    eik_solver.cycFIM();
    eik_solver.clean_list();

    eik_solver.stats.log_stats(user_globals::tm_manager->time, false);
    eik_solver.nodeData.log_stats(user_globals::tm_manager->time, false);
    do_output_eikonal = true;
  }
}

void Eikonal::solve_RD()
{
  // if requested, we checkpoint the current state
  checkpointing();

  // activation checking
  const double time      = user_globals::tm_manager->time,
               time_step = user_globals::tm_manager->time_step;

  lat.check_acts(time);
  lat.check_quiescence(time, time_step);

  clamp_Vm();

  *parab_solver.old_vm = *parab_solver.Vmv;  // needed in step D of DREAM

  // compute ionics update
  ion.compute_step();

  // Compute the I diff current
  *eik_solver.Idiff *= 0;
  eik_solver.compute_diffusion_current(time, *parab_solver.Vmv);
  parab_solver.Vmv->add_scaled(*eik_solver.Idiff, time_step);

  if (eik_tech == REp) {
    parab_solver.solve(*phie_dummy);
  }

  switch (eik_tech) {
    case DREAM: eik_solver.update_repolarization_times_from_rd(*parab_solver.Vmv, *parab_solver.old_vm, time); break;
    default: break;
  }

  clamp_Vm();
}

void eikonal_solver::init()
{
  if (param_globals::output_level > 1) log_msg(0, 0, 0, "\n    *** Initializing Eikonal Solver ***\n");

  if (get_size() > 1) {
    log_msg(NULL, 5, 0, "The DREAM does not run in parallel. Aborting!");
    EXIT(EXIT_FAILURE);
  }

  stats.init_logger("eik_stats.dat");

  if (param_globals::dream.output.debugNode >= 0) {
    nodeData.idX = param_globals::dream.output.debugNode;
    char buf[256];
    snprintf(buf, sizeof buf, "node_%d.dat", nodeData.idX);
    nodeData.init_logger(buf);
  }

  // currently only used for output
  const sf_mesh& mesh = get_mesh(intra_elec_msh);
  const int      dpn  = 1;
  SF::init_vector(&Idiff, mesh, dpn, sf_vec::algebraic);
  SF::init_vector(&AT, mesh, dpn, sf_vec::algebraic);
  SF::init_vector(&RT, mesh, dpn, sf_vec::algebraic);

  num_pts  = mesh.l_numpts;

  e2n_con = mesh.con;     // Connectivity vector with nodes that belong to each element
                          // The 4 nodes from index 4j to 4j+3 belong to the same element for 0<=j<num elements
  elem_start = mesh.dsp;  // For the i_th element elem_start[j] stores its starting position in the "connect" vector

  cnt_from_dsp(elem_start, e2n_cnt);
  transpose_connectivity(e2n_cnt, e2n_con, n2e_cnt, n2e_con);
  dsp_from_cnt(n2e_cnt, n2e_dsp);

  x_coord.resize(num_pts);
  y_coord.resize(num_pts);
  z_coord.resize(num_pts);
  fiber_node.resize(num_pts * 3, 0);
  sheet_node.resize(num_pts * 3, 0);

  for (int i = 0; i < ((mesh.xyz.size() + 1) / 3); i++) {
    x_coord[i] = mesh.xyz[3 * i + 0];
    y_coord[i] = mesh.xyz[3 * i + 1];
    z_coord[i] = mesh.xyz[3 * i + 2];
  }

  diff_cur.resize(num_pts);

  D.set_size(3, 3);
  iD.set_size(3, 3);

  List.resize(num_pts, 0);           // List of active nodes
  T_A.resize(num_pts, inf);          // Activation times
  D_I.resize(num_pts, 0);            // Diastolic Intervals
  T_R.resize(num_pts, -400);         // Recovery times
  TA_old.resize(num_pts, inf);       // Activation Time is the previous cycle
  num_changes.resize(num_pts, 0);    // Number of Times entered in the List per current LAT
  nReadded2List.resize(num_pts, 0);  // Number of Times entered in the List per current LAT
  stim_status.resize(num_pts, 0);    // Status of nodes

  CVnodes.resize(num_pts, 0);
  p_cvrest.resize(num_pts);
  k_cvrest.resize(num_pts);
  j_cvrest.resize(num_pts);
  l_cvrest.resize(num_pts);
  cv_ar_t.resize(num_pts, -1);
  cv_ar_n.resize(num_pts, -1);

  StimulusPoints.resize(num_pts * 10, -1);  // Nodes that have initial stimulus
  StimulusTimes.resize(num_pts * 10, -1);   // Times when the initial nodes are stimulated

  switch (mesh.type[0]) {
    case 0:
      numPtsxElem = 4;
      break;

    case 6:
      numPtsxElem = 3;
      break;

    case 7:
      numPtsxElem = 2;
      break;

    default:
      log_msg(0, 5, 0, "Error: Type of element is not compatible with this version of the eikonal model. Use tetrahedra or triangles");
      EXIT(EXIT_FAILURE);
      break;
  }

  if (strlen(param_globals::start_statef) > 0) load_state_file();

  create_node_to_node_graph();

  translate_stim_to_eikonal();

  init_diffusion_current();
}

void eikonal_solver::bubble_sort(SF::vector<SF_real>& Array, SF::vector<mesh_int_t>& IndexArray, int& length)
{
  SF_real    Temp_a;
  mesh_int_t Temp_i;
  for (int i = 0; i < length - 1; i++) {
    for (int j = 0; j < length - i - 1; j++) {
      if (Array[j] > Array[j + 1]) {
        Temp_a            = Array[j];
        Array[j]          = Array[j + 1];
        Array[j + 1]      = Temp_a;
        Temp_i            = IndexArray[j];
        IndexArray[j]     = IndexArray[j + 1];
        IndexArray[j + 1] = Temp_i;
      }
    }
  }
}

bool eikonal_solver::determine_model_to_run(double& time)
{
  bool act_is_in_safety_window, empty_list_and_illegal_stimulus, empty_stimulus, stimulus_is_in_safety_window;

  act_is_in_safety_window         = (actMIN > param_globals::dream.tau_s) && (actMIN - time > param_globals::dream.tau_s);
  empty_stimulus                  = (StimulusPoints[Index_currStim] == -1);
  stimulus_is_in_safety_window    = (StimulusTimes[Index_currStim] - time > param_globals::dream.tau_s);
  empty_list_and_illegal_stimulus = (sum(List) == 0) && (empty_stimulus || stimulus_is_in_safety_window);

  return act_is_in_safety_window || empty_list_and_illegal_stimulus;
}

void eikonal_solver::set_initial_cv()
{
  const sf_mesh&         mesh = get_mesh(intra_elec_msh);
  SF::vector<mesh_int_t> vertices;
  for (int num_g = 0; num_g < param_globals::num_gregions; num_g++) {
    for (int num_id = 0; num_id <= param_globals::gregion[num_g].num_IDs; num_id++) {
      indices_from_region_tag(vertices, mesh, param_globals::gregion[num_g].ID[num_id]);
      for (int i = 0; i < vertices.size(); i++) {
        CVnodes[vertices[i]] = param_globals::gregion[num_g].dream.vel_l;
        cv_ar_t[vertices[i]] = param_globals::gregion[num_g].dream.vel_l / param_globals::gregion[num_g].dream.vel_t;
        if (twoFib) {
          cv_ar_n[vertices[i]] = param_globals::gregion[num_g].dream.vel_l / param_globals::gregion[num_g].dream.vel_n;
        }
      }
    }
  }
  CVnodes_mod = CVnodes;

  for (int num_i = 0; num_i < param_globals::num_imp_regions; num_i++) {
    for (int num_id = 0; num_id <= param_globals::imp_region[num_i].num_IDs; num_id++) {
      indices_from_region_tag(vertices, mesh, param_globals::imp_region[num_i].ID[num_id]);
      for (int i = 0; i < vertices.size(); i++) {
        p_cvrest[vertices[i]] = param_globals::imp_region[num_i].dream.CVrest.rho;
        k_cvrest[vertices[i]] = param_globals::imp_region[num_i].dream.CVrest.kappa;
        j_cvrest[vertices[i]] = param_globals::imp_region[num_i].dream.CVrest.theta;
        l_cvrest[vertices[i]] = param_globals::imp_region[num_i].dream.CVrest.psi;
      }
    }
  }

  if (strlen(param_globals::cv_scale_vec)) {
    mesh_t      mt   = intra_elec_msh;
    const char* file = param_globals::cv_scale_vec;

    size_t num_file_entries = SF::root_count_ascii_lines(file, mesh.comm);

    if (num_file_entries != mesh.g_numpts)
      log_msg(0, 2, 0, "%s warning: number of %s conductivity scaling entries does not match number of points!,%i,%i",
              __func__, get_mesh_type_name(mt), num_file_entries, mesh.g_numpts);

    // set up parallel point vector and read data
    sf_vec* escale;
    SF::init_vector(&escale, mesh, 1, sf_vec::nodal);
    escale->read_ascii(file);

    if (get_size() > 1) {
      // set up point vector permutation and permute
      if (get_permutation(mt, PETSC_TO_CANONICAL, 1) == NULL) {
        register_permutation(mt, PETSC_TO_CANONICAL, 1);
      }
      SF::scattering& sc = *get_permutation(mt, PETSC_TO_CANONICAL, 1);
      sc(*escale, false);
    }

    // copy data into SF::vector
    SF::vector<double> nd_scale;
    SF_real*           p = escale->ptr();
    nd_scale.assign(p, p + escale->lsize());
    escale->release_ptr(p);

    for (int num_g = 0; num_g < param_globals::num_gregions; num_g++) {
      for (int num_id = 0; num_id < param_globals::gregion[num_g].num_IDs; num_id++) {
        indices_from_region_tag(vertices, mesh, param_globals::gregion[num_g].ID[num_id]);
        for (int i = 0; i < vertices.size(); i++) {
          CVnodes[vertices[i]] *= nd_scale[vertices[i]];
        }
      }
    }
    CVnodes_mod = CVnodes;
  }
}  // Close set initial CV

void eikonal_solver::update_Ta_in_active_list()
{
  bool First_in_List = 1;

  for (size_t j = 0; j < List.size(); j++) {
    if (T_A[j] < 0) continue;
    if (List[j] == 1 && First_in_List) {
      actMIN        = T_A[j];
      actMAX        = T_A[j];
      First_in_List = 0;
      continue;
    }
    if (List[j] == 1) {
      if (T_A[j] < actMIN) actMIN = T_A[j];
      if (T_A[j] > actMAX) actMAX = T_A[j];
    }
  }
}

void eikonal_solver::clean_list()
{
  for (size_t j = 0; j < List.size(); j++) {
    if (T_A[j] < user_globals::tm_manager->time) {
      TA_old[j]        = T_A[j];
      stim_status[j]   = 0;
      nReadded2List[j] = 0;

      if (List[j] == 1) List[j] = 0;
    }
  }
}

SF_real eikonal_solver::compute_H(mesh_int_t& indX, SF_real& T_A_indX)
{
  float delay = 0;
  if (stim_status[indX] == 2) delay = 5;

  if (T_A_indX == inf || (T_R[indX] + delay) == -400 || abs(T_A_indX - (TA_old[indX])) < 2) {
    return CVnodes[indX];
  } else if ((T_R[indX]) > T_A_indX) {
    return CVnodes[indX];
  } else {
    float DI_indX = T_A_indX - (T_R[indX] + delay);
    D_I[indX]     = DI_indX;

    float Factor_DI;
    float denom = log(p_cvrest[indX]) / l_cvrest[indX];
    Factor_DI   = 1 - p_cvrest[indX] * exp(-(DI_indX + k_cvrest[indX]) * denom);

    if (Factor_DI > 0) {
      if (DI_indX < j_cvrest[indX]) {
        return -CVnodes[indX] * Factor_DI;
      } else {
        return CVnodes[indX] * Factor_DI;
      }
    } else {
      return 0;
    }
  }
}

SF_real eikonal_solver::compute_coherence(mesh_int_t& indX)
{
  SF_real p       = T_A[indX];
  int     iter_ch = 0;
  SF_real q       = compute_F(indX, p);
  while (abs(p - q) >= param_globals::dream.fim.tol && iter_ch < param_globals::dream.fim.max_coh) {
    p = q;
    q = compute_F(indX, p);
    iter_ch++;
  }
  if (compute_H(indX, q) < 0) return T_A[indX];
  return q;
}

void eikonal_solver::compute_bc()
{
  bool EmpList = sum(List) == 0;

  if (EmpList) {
    for (size_t j = Index_currStim; j < StimulusTimes.size(); j++) {
      if ((StimulusPoints[j] == -1) || (StimulusTimes[j] > StimulusTimes[Index_currStim])) {
        Index_currStim = j;
        break;
      }

      mesh_int_t indNode = StimulusPoints[j];
      SF_real    TimeSt  = StimulusTimes[j];

      bool cond2add = add_node_neighbor_to_list(T_R[indNode], T_A[indNode], TimeSt);
      if (cond2add && compute_H(indNode, TimeSt) > 0) {
        if (T_A[indNode] < TimeSt) {
          CVnodes_mod[indNode] = -CVnodes_mod[indNode];
        }
        T_A[indNode]         = TimeSt;
        stim_status[indNode] = 1;
        stats.bc_status      = true;

        if (param_globals::dream.output.debugNode == indNode) {
          nodeData.T_A = TimeSt;
          nodeData.update_status(node_stats::in, node_stats::stim);
        }

        for (int ii = n2n_dsp[indNode]; ii < n2n_dsp[indNode + 1]; ii++) {
          mesh_int_t indXNB = n2n_connect[ii];
          if (List[indXNB] == 0) {
            List[indXNB] = 1;
            if (param_globals::dream.output.debugNode == indXNB) {
              nodeData.update_status(node_stats::in, node_stats::nbn);
              nodeData.idXNB   = indNode;
              nodeData.nbn_T_A = TimeSt;
            }
          }
        }
      }

      if (cond2add && compute_H(indNode, TimeSt) <= 0) {
        stim_status[indNode] = 2;
      }
    }

    // After NB of initial points are assigned remove initial points from the list if they were added in the previous loop.
    for (size_t j = 0; j < StimulusPoints.size(); j++) {
      if ((StimulusPoints[j] == -1) && (StimulusTimes[j] > StimulusTimes[Index_currStim])) continue;
      if (List[StimulusPoints[j]] == 1) {
        List[StimulusPoints[j]] = 0;
        if (param_globals::dream.output.debugNode == StimulusPoints[j]) nodeData.update_status(node_stats::out, node_stats::stim);
      }
    }

    if (StimulusTimes[Index_currStim] != -1) actMAX = StimulusTimes[Index_currStim];
  } else {
    for (size_t i = Index_currStim; i < StimulusTimes.size(); i++) {
      mesh_int_t indNode = StimulusPoints[i];
      SF_real    TimeSt  = StimulusTimes[i];

      if (indNode == -1) continue;

      if (TimeSt <= actMAX) {
        bool cond2add = add_node_neighbor_to_list(T_R[indNode], T_A[indNode], TimeSt);
        if (cond2add && compute_H(indNode, TimeSt) > 0) {
          T_A[indNode]         = TimeSt;
          stim_status[indNode] = 1;
          stats.bc_status      = true;

          if (param_globals::dream.output.debugNode == indNode) {
            nodeData.T_A = TimeSt;
            nodeData.update_status(node_stats::in, node_stats::stim);
          }

          for (int ii = n2n_dsp[indNode]; ii < n2n_dsp[indNode + 1]; ii++) {
            mesh_int_t indXNB = n2n_connect[ii];
            if (List[indXNB] == 0) {
              List[indXNB] = 1;
              if (param_globals::dream.output.debugNode == indXNB) {
                nodeData.update_status(node_stats::in, node_stats::nbn);
                nodeData.idXNB   = indNode;
                nodeData.nbn_T_A = TimeSt;
              }
            }
          }

          for (size_t j = 0; j < StimulusPoints.size(); j++) {
            if ((StimulusPoints[j] == -1) && (StimulusTimes[j] > StimulusTimes[Index_currStim])) continue;
            if (List[StimulusPoints[j]] == 1) {
              List[StimulusPoints[j]] = 0;
              if (param_globals::dream.output.debugNode == StimulusPoints[j]) nodeData.update_status(node_stats::out, node_stats::stim);
            }
          }
        }
        if (cond2add && compute_H(indNode, TimeSt) <= 0) {
          stim_status[indNode] = 2;
        }

      } else {
        Index_currStim = i;

        break;
      }
    }
  }

}  // end compute stimulus

void eikonal_solver::FIM()
{
  int    niter = 0;
  double t1, t0;
  get_time(t0);

  for (size_t i = Index_currStim; i < StimulusTimes.size(); i++) {
    mesh_int_t indNode = StimulusPoints[i];
    SF_real    TimeSt  = StimulusTimes[i];
    if (StimulusPoints[i] == -1) break;

    T_A[indNode] = TimeSt;

    for (size_t j = n2n_dsp[indNode]; j < n2n_dsp[indNode + 1]; j++) {
      if (T_A[n2n_connect[j]] == inf) {
        List[n2n_connect[j]] = 1;
      }
    }
  }

  do {
    for (mesh_int_t indX = 0; indX < List.size(); indX++) {
      if (List[indX] == 0) continue;

      SF_real p = T_A[indX];
      SF_real q = compute_F(indX, p);
      T_A[indX] = q;

      if (abs(p - q) < param_globals::dream.fim.tol) {
        for (int ii = n2n_dsp[indX]; ii < n2n_dsp[indX + 1]; ii++) {
          mesh_int_t indXNB = n2n_connect[ii];
          if (List[indXNB] == 1) continue;
          p = T_A[indXNB];
          q = compute_F(indXNB, p);
          if (p > q) {
            T_A[indXNB]  = q;
            List[indXNB] = 1;
          }
        }
        List[indX] = 0;
      }
    }

    update_Ta_in_active_list();
    niter++;

    if (niter > param_globals::dream.fim.max_iter) break;

  } while (sum(List) > 0);

  // copy for igb output
  double* atc = AT->ptr();
  for (mesh_int_t i = 0; i < List.size(); i++) {
    if (T_A[i] == inf) {
      // for better visualization in meshalyzer
      atc[i] = -1;
    } else {
      atc[i] = T_A[i];
    }
  }
  AT->release_ptr(atc);

  // treat solver statistics
  auto dur = timing(t1, t0);
  stats.slvtime_A += dur;
  stats.update_iter(niter);
  stats.minAT      = actMIN;
  stats.maxAT      = actMAX;
  stats.activeList = sum(List);
  stats.update_cli(user_globals::tm_manager->time, false);
}

void eikonal_solver::cycFIM()
{
  int    niter = 0;
  double t1, t0;
  get_time(t0);

  if (sum(List) == 0) {
    compute_bc();
  }

  time2stop_eikonal = param_globals::dream.tau_inc;
  float maxadvance  = user_globals::tm_manager->time + param_globals::dream.tau_s + param_globals::dream.tau_inc + param_globals::dream.tau_max;

  do {
    if (StimulusTimes[Index_currStim] <= maxadvance) {
      compute_bc();
    }

    for (mesh_int_t indX = 0; indX < List.size(); indX++) {
      SF_real p = T_A[indX];
      SF_real q;

      if (List[indX] == 0) continue;

      q = compute_coherence(indX);

      T_A[indX]         = q;
      num_changes[indX] = num_changes[indX] + 1;

      if (q > maxadvance) continue;

      if (abs(p - q) < param_globals::dream.fim.tol || (num_changes[indX] > param_globals::dream.fim.max_iter) || (q == inf || p == inf)) {
        for (int ii = n2n_dsp[indX]; ii < n2n_dsp[indX + 1]; ii++) {
          mesh_int_t indXNB = n2n_connect[ii];

          if (List[indXNB] == 1) {
            continue;
          }

          SF_real pNB = T_A[indXNB];
          SF_real qNB;

          qNB = compute_coherence(indXNB);

          bool node_is_valid = add_node_neighbor_to_list(T_R[indXNB], pNB, qNB) && qNB != inf && qNB > user_globals::tm_manager->time;
          // This second condition is there to be able to add a node to the list when a new valid activation time
          // is found but a reentry would be blocked by the L2 parameter. Normally this condition does not make or brake the
          // simulation but would leave individual nodes inactivated, which is not ideal.
          bool ignore_L2_if_valid = node_is_valid && qNB > pNB && nReadded2List[indXNB] >= param_globals::dream.fim.max_addpt;

          if (node_is_valid && (nReadded2List[indXNB] < param_globals::dream.fim.max_addpt) || ignore_L2_if_valid) {
            if (qNB > pNB) {
              nReadded2List[indXNB] = 0;
            }
            T_A[indXNB] = qNB;
            nReadded2List[indXNB]++;
            num_changes[indXNB] = 0;
            List[indXNB]        = 1;
            if (param_globals::dream.output.debugNode == indXNB) {
              nodeData.update_status(node_stats::in, node_stats::nbn);
              nodeData.idXNB   = indX;
              nodeData.nbn_T_A = q;
            }
          }
        }

        List[indX] = 0;
        if (param_globals::dream.output.debugNode == indX) {
          nodeData.update_status(node_stats::out, node_stats::conv);
        }
      }
    }

    SF_real actMIN_old = actMIN;
    SF_real progress_time;

    update_Ta_in_active_list();

    if (actMIN > actMIN_old) {
      progress_time = actMIN - actMIN_old;
    } else {
      progress_time = 0;
    }

    time2stop_eikonal -= progress_time;
    niter++;

  } while (time2stop_eikonal > 0 && sum(List) > 0);

  if (sum(List) == 0) {
    if (StimulusTimes[Index_currStim] != -1) actMIN = StimulusTimes[Index_currStim];
  }

  // copy for igb output
  double* atc = AT->ptr();
  double* rpt = RT->ptr();
  for (mesh_int_t i = 0; i < List.size(); i++) {
    if (T_A[i] == inf) {
      // for better visualization in meshalyzer
      atc[i] = -1;
    } else {
      atc[i] = T_A[i];
    }
    rpt[i] = T_R[i];
  }
  AT->release_ptr(atc);
  RT->release_ptr(rpt);

  // treat solver statistics
  auto dur = timing(t1, t0);
  stats.slvtime_A += dur;
  stats.update_iter(niter);
  stats.minAT      = actMIN;
  stats.maxAT      = actMAX;
  stats.activeList = sum(List);
  stats.update_cli(user_globals::tm_manager->time, false);

  // treat node stats
  if (param_globals::dream.output.debugNode >= 0) {
    nodeData.T_A = T_A[nodeData.idX];
    nodeData.T_R = T_R[nodeData.idX];
    nodeData.D_I = D_I[nodeData.idX];
  }

}  // close iterate list

SF_real eikonal_solver::compute_F(mesh_int_t& indX, SF_real& T_A_indX)
{
  const double time = user_globals::tm_manager->time;
  TA                = inf;
  SF::vector<mesh_int_t> Winner_Element(3, -1);
  CVnodes_mod[indX] = compute_H(indX, T_A_indX);

  CV_l = abs(CVnodes_mod[indX]);

  if (CV_l == 0) {
    return TA;
  }

  SF::Point fib, she;

  if (twoFib) {
    fib.x = fiber_node[3 * indX + 0];
    fib.y = fiber_node[3 * indX + 1];
    fib.z = fiber_node[3 * indX + 2];

    she.x = sheet_node[3 * indX + 0];
    she.y = sheet_node[3 * indX + 1];
    she.z = sheet_node[3 * indX + 2];

    SF_real   ani_ratio_t2 = cv_ar_t[indX] * cv_ar_t[indX];
    SF_real   ani_ratio_n2 = cv_ar_n[indX] * cv_ar_n[indX];
    SF::Point n            = cross(fib, she);
    SF::outer_prod(fib, fib, 1, D[0], false);
    SF::outer_prod(she, she, 1 / ani_ratio_t2, D[0], true);
    SF::outer_prod(n, n, 1 / ani_ratio_n2, D[0], true);

    iD = invert_3x3(D);

  } else {
    fib.x = fiber_node[3 * indX + 0];
    fib.y = fiber_node[3 * indX + 1];
    fib.z = fiber_node[3 * indX + 2];

    SF_real ani_ratio2 = cv_ar_t[indX] * cv_ar_t[indX];
    SF::outer_prod(fib, fib, (ani_ratio2 - 1) / ani_ratio2, D[0]);
    D[0][0] += (1 / ani_ratio2), D[1][1] += (1 / ani_ratio2), D[2][2] += (1 / ani_ratio2);
    iD = invert_3x3(D);
  }

  if (numPtsxElem == 4) {
    // Solve for Tetrahedron
    for (int e_i = n2e_dsp[indX]; e_i < n2e_dsp[indX + 1]; e_i++) {
      int        Elem_i = n2e_con[e_i];
      mesh_int_t indEle = elem_start[Elem_i];

      std::vector<mesh_int_t> Elem_cur = {e2n_con[indEle], e2n_con[indEle + 1], e2n_con[indEle + 2], e2n_con[indEle + 3]};
      std::remove(Elem_cur.begin(), Elem_cur.end(), indX);  // intentionally used to move indX to the end of Elem_cur but not to be erased

      condign1 = T_A[Elem_cur[0]] == inf || T_A[Elem_cur[1]] == inf || T_A[Elem_cur[2]] == inf || (T_A[Elem_cur[0]] < time) || (T_A[Elem_cur[1]] < time) || (T_A[Elem_cur[2]] < time);
      condign2 = (T_A[Elem_cur[0]] < T_R[Elem_cur[0]]) || (T_A[Elem_cur[1]] < T_R[Elem_cur[1]]) || (T_A[Elem_cur[2]] < T_R[Elem_cur[2]]);
      condign3 = CV_l == 0;
      condign4 = stim_status[indX] == 2 && (stim_status[Elem_cur[0]] == 1 || stim_status[Elem_cur[1]] == 1 || stim_status[Elem_cur[2]] == 1);

      if (!(condign1 || condign2 || condign3 || condign4)) {
        SF_real temp3 = solve_tet(indX, Elem_cur);

        if (temp3 < TA && temp3 > T_R[indX] && compute_H(indX, temp3) > 0) {
          Winner_Element[0] = Elem_cur[0];
          Winner_Element[1] = Elem_cur[1];
          Winner_Element[2] = Elem_cur[2];
          TA                = temp3;
        }
      }

      std::vector<mesh_int_t> Triang = {-1, -1};
      // Solve triangles in Tetrahedra
      for (int cases = 0; cases < 3; cases++) {
        switch (cases) {
          case 0:
            Triang[0] = Elem_cur[0];
            Triang[1] = Elem_cur[1];
            break;

          case 1:
            Triang[0] = Elem_cur[0];
            Triang[1] = Elem_cur[2];
            break;

          case 2:
            Triang[0] = Elem_cur[1];
            Triang[1] = Elem_cur[2];
            break;

          default:
            break;
        }

        condign1 = T_A[Triang[0]] == inf || T_A[Triang[1]] == inf || (T_A[Triang[0]] < time) || (T_A[Triang[1]] < time);
        condign2 = (T_A[Triang[0]] < T_R[Triang[0]]) || (T_A[Triang[1]] < T_R[Triang[1]]);  //||abs(T_A[Triang[0]] - TA_old[Triang[0]])<2||abs(T_A[Triang[1]] - TA_old[Triang[1]])<2;
        condign3 = CV_l == 0;
        condign4 = stim_status[indX] == 2 && (stim_status[Triang[0]] == 1 || stim_status[Triang[1]] == 1);

        if (!(condign1 || condign2 || condign3 || condign4)) {
          SF_real temp2 = solve_tri(indX, Triang);

          if (temp2 < TA && temp2 > T_R[indX] && compute_H(indX, temp2) > 0) {
            Winner_Element[0] = Triang[0];
            Winner_Element[1] = Triang[1];
            Winner_Element[2] = -1;
            TA                = temp2;
          }
        }
      }
    }
  }  // numPtsxElem == 4

  if (numPtsxElem == 3) {
    for (int e_i = n2e_dsp[indX]; e_i < n2e_dsp[indX + 1]; e_i++) {
      int        Elem_i = n2e_con[e_i];
      mesh_int_t indEle = elem_start[Elem_i];

      std::vector<mesh_int_t> Elem_cur = {e2n_con[indEle], e2n_con[indEle + 1], e2n_con[indEle + 2]};
      std::remove(Elem_cur.begin(), Elem_cur.end(), indX);  // intentionally used to move indX to the end of Elem_cur but not to be erased

      condign1 = T_A[Elem_cur[0]] == inf || T_A[Elem_cur[1]] == inf || (T_A[Elem_cur[0]] < time) || (T_A[Elem_cur[1]] < time);  //||abs(T_A[Elem_cur[0]] - TA_old[Elem_cur[0]])<2||abs(T_A[Elem_cur[1]] - TA_old[Elem_cur[1]])<2;
      condign2 = (T_A[Elem_cur[0]] < T_R[Elem_cur[0]]) || (T_A[Elem_cur[1]] < T_R[Elem_cur[1]]);
      condign3 = CV_l == 0;
      condign4 = stim_status[indX] == 2 && (stim_status[Elem_cur[0]] == 1 || stim_status[Elem_cur[1]] == 1);

      if (condign1 || condign2 || condign3 || condign4) continue;

      SF_real temp2 = solve_tri(indX, Elem_cur);

      if (TA > temp2 && temp2 > T_R[indX] && compute_H(indX, temp2) > 0) {
        TA = temp2;
      }
    }
  }  // numPtsxElem == 3

  // Solve for NB Nodes
  for (int ii = n2n_dsp[indX]; ii < n2n_dsp[indX + 1]; ii++) {
    mesh_int_t indXNB = n2n_connect[ii];
    condign4          = stim_status[indX] == 2 && (stim_status[indXNB] == 1);

    if (!(CV_l == 0 || T_A[indXNB] < T_R[indXNB] || T_A[indXNB] == inf || condign4 || T_A[indXNB] < time)) {
      SF_real temp1 = solve_edges(indX, indXNB);

      if (temp1 < TA && temp1 > T_R[indX] && compute_H(indX, temp1) > 0) {
        Winner_Element[0] = indXNB;
        Winner_Element[1] = -1;
        Winner_Element[2] = -1;
        TA                = temp1;
      }
    }
  }

  if (TA == inf) {
    Winner_Element[0] = -1;
    Winner_Element[1] = -1;
    Winner_Element[2] = -1;
    TA                = T_A[indX];
  }

  std::sort(Winner_Element.begin(), Winner_Element.end());
  if (compute_H(indX, TA) == 0) TA = inf;

  return TA;

}  // Close Function_F

SF_real eikonal_solver::solve_tet(mesh_int_t& indx3, std::vector<mesh_int_t>& Elem_cur)
{
  SF_real TA_th = inf;

  SF::vector<SF_real> X(3), Y(3), Z(3), W(3);
  SF::vector<SF_real> XY(3), XZ(3), XW(3);

  int indy3 = Elem_cur[0], indz3 = Elem_cur[1], indw3 = Elem_cur[2];

  X[0] = x_coord[indx3];
  X[1] = y_coord[indx3];
  X[2] = z_coord[indx3];

  Y[0] = x_coord[indy3];
  Y[1] = y_coord[indy3];
  Y[2] = z_coord[indy3];

  Z[0] = x_coord[indz3];
  Z[1] = y_coord[indz3];
  Z[2] = z_coord[indz3];

  W[0] = x_coord[indw3];
  W[1] = y_coord[indw3];
  W[2] = z_coord[indw3];

  XY[0] = Y[0] - X[0];
  XZ[0] = Z[0] - X[0];
  XW[0] = W[0] - X[0];
  XY[1] = Y[1] - X[1];
  XZ[1] = Z[1] - X[1];
  XW[1] = W[1] - X[1];
  XY[2] = Y[2] - X[2];
  XZ[2] = Z[2] - X[2];
  XW[2] = W[2] - X[2];

  c1 = mult_VtMV(XY, iD, XY);
  c2 = mult_VtMV(XZ, iD, XZ);
  c3 = mult_VtMV(XW, iD, XW);
  c4 = 2 * (mult_VtMV(XY, iD, XZ));  //+mult_VtMV(XZ,iD,XY);
  c5 = 2 * mult_VtMV(XY, iD, XW);    //+mult_VtMV(XW,iD,XY);
  c6 = 2 * mult_VtMV(XW, iD, XZ);    //+mult_VtMV(XZ,iD,XW);

  a1          = 4 * c2 * c3 - 2 * c2 * c5 - 2 * c3 * c4 + c4 * c6 + c5 * c6 - c6 * c6;
  s1          = 4 * c1 * c2 + 4 * c1 * c3 + 4 * c2 * c3 - 4 * c1 * c6 - 4 * c2 * c5 - 4 * c3 * c4 + 2 * c4 * c5 + 2 * c4 * c6 + 2 * c5 * c6 - c4 * c4 - c5 * c5 - c6 * c6;
  SF_real s_1 = 1 / s1;
  s2          = c1 * c6 * c6 + c2 * c5 * c5 + c3 * c4 * c4 - 4 * c1 * c2 * c3 - c4 * c5 * c6;
  s3          = 4 * c1 * c6 - 4 * c1 * c3 - 4 * c2 * c3 - 4 * c1 * c2 + 4 * c2 * c5 + 4 * c3 * c4 - 2 * c4 * c5 - 2 * c4 * c6 - 2 * c5 * c6 + c4 * c4 + c5 * c5 + c6 * c6 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indy3] * c2 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indy3] * c3 - 4 * CV_l * CV_l * T_A[indy3] * T_A[indy3] * c6 + 4 * CV_l * CV_l * T_A[indz3] * T_A[indz3] * c1 + 4 * CV_l * CV_l * T_A[indz3] * T_A[indz3] * c3 - 4 * CV_l * CV_l * T_A[indz3] * T_A[indz3] * c5 + 4 * CV_l * CV_l * T_A[indw3] * T_A[indw3] * c1 + 4 * CV_l * CV_l * T_A[indw3] * T_A[indw3] * c2 - 4 * CV_l * CV_l * T_A[indw3] * T_A[indw3] * c4 - 8 * CV_l * CV_l * T_A[indy3] * T_A[indz3] * c3 - 4 * CV_l * CV_l * T_A[indy3] * T_A[indz3] * c4 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indz3] * c5 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indz3] * c6 - 8 * CV_l * CV_l * T_A[indy3] * T_A[indw3] * c2 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indw3] * c4 - 4 * CV_l * CV_l * T_A[indy3] * T_A[indw3] * c5 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indw3] * c6 - 8 * CV_l * CV_l * T_A[indz3] * T_A[indw3] * c1 + 4 * CV_l * CV_l * T_A[indz3] * T_A[indw3] * c4 + 4 * CV_l * CV_l * T_A[indz3] * T_A[indw3] * c5 - 4 * CV_l * CV_l * T_A[indz3] * T_A[indw3] * c6;
  a2          = 2 * CV_l * T_A[indy3] * c2 + 2 * CV_l * T_A[indy3] * c3 - 2 * CV_l * T_A[indy3] * c6 - 2 * CV_l * T_A[indz3] * c3 - CV_l * T_A[indz3] * c4 + CV_l * T_A[indz3] * c5 + CV_l * T_A[indz3] * c6 - 2 * CV_l * T_A[indw3] * c2 + CV_l * T_A[indw3] * c4 - CV_l * T_A[indw3] * c5 + CV_l * T_A[indw3] * c6;
  s4          = 4 * c1 * c2 + 4 * c1 * c3 + 4 * c2 * c3 - 4 * c1 * c6 - 4 * c2 * c5 - 4 * c3 * c4 + 2 * c4 * c5 + 2 * c4 * c6 + 2 * c5 * c6 - c4 * c4 - c5 * c5 - c6 * c6;
  b1          = 4 * c1 * c3 - 2 * c1 * c6 - 2 * c3 * c4 + c4 * c5 + c5 * c6 - c5 * c5;
  b2          = CV_l * T_A[indy3] * c5 - CV_l * T_A[indy3] * c4 - 2 * CV_l * T_A[indy3] * c3 + CV_l * T_A[indy3] * c6 + 2 * CV_l * T_A[indz3] * c1 + 2 * CV_l * T_A[indz3] * c3 - 2 * CV_l * T_A[indz3] * c5 - 2 * CV_l * T_A[indw3] * c1 + CV_l * T_A[indw3] * c4 + CV_l * T_A[indw3] * c5 - CV_l * T_A[indw3] * c6;

  SF_real s2_s3 = s2 / s3;
  SF_real LTas  = a1 * s_1;
  SF_real LTbs  = b1 * s_1;

  SF_real RTa = (2 * sqrt(s2_s3) * (a2)) / (s4);
  SF_real RTb = (2 * sqrt(s2_s3) * (b2)) / (s4);

  SF_real tol = 1e-24;

  if (abs(s1) > tol && abs(s3) > tol && abs(s4) > tol && (s2 / s3) > 0) {
    a = LTas - RTa;
    b = LTbs - RTb;

    if ((a > 0) && (a < 1) && (b > 0) && (b < 1) && ((1 - a - b) > 0) && ((1 - a - b) < 1)) {
      TA_th = std::min(TA_th, compute_AT_at_node_in_tet(a, b, indx3, indy3, indz3, indw3));
    }

    a = LTas + RTa;
    b = LTbs + RTb;

    if ((a > 0) && (a < 1) && (b > 0) && (b < 1) && ((1 - a - b) > 0) && ((1 - a - b) < 1)) {
      TA_th = std::min(TA_th, compute_AT_at_node_in_tet(a, b, indx3, indy3, indz3, indw3));
    }
  }

  return TA_th;
}

SF_real eikonal_solver::solve_tri(mesh_int_t& indx, std::vector<mesh_int_t>& Elem_cur)
{
  SF_real TA_tr = inf;

  SF::vector<SF_real> X(3), Y(3), Z(3);
  SF::vector<SF_real> XY(3), XZ(3);

  int indy = Elem_cur[0], indz = Elem_cur[1];

  X[0] = x_coord[indx];
  X[1] = y_coord[indx];
  X[2] = z_coord[indx];

  Y[0] = x_coord[indy];
  Y[1] = y_coord[indy];
  Y[2] = z_coord[indy];

  Z[0] = x_coord[indz];
  Z[1] = y_coord[indz];
  Z[2] = z_coord[indz];

  XY[0] = Y[0] - X[0];
  XZ[0] = Z[0] - X[0];
  XY[1] = Y[1] - X[1];
  XZ[1] = Z[1] - X[1];
  XY[2] = Y[2] - X[2];
  XZ[2] = Z[2] - X[2];

  C = mult_VtMV(XZ, iD, XZ);
  A = mult_VtMV(XY, iD, XY) - 2 * mult_VtMV(XZ, iD, XY) + C;
  B = 2 * mult_VtMV(XY, iD, XZ) - 2 * C;

  if (abs(T_A[indy] - T_A[indz]) < 1e-5) {
    p1 = -B / (2 * A);
    if (p1 > 0 && p1 < 1) {
      TA_tr = std::min(compute_AT_at_node(p1, indx, indy, indz), TA_tr);
    }
  }

  lam = ((T_A[indz] - T_A[indy]) * (T_A[indz] - T_A[indy])) * 4 * (CV_l * CV_l);
  ae  = (4 * (A * A)) - lam * A;
  be  = 4 * A * B - lam * B;
  ce  = (B * B) - lam * C;
  det = (be * be) - (4 * ae * ce);

  if (det >= 0) {
    p1 = (-be + sqrt(det)) / (2 * ae);

    if (p1 > 0 && p1 < 1) {
      TA_tr = std::min(compute_AT_at_node(p1, indx, indy, indz), TA_tr);
    }

    if (det > 0) {
      p2 = (-be - sqrt(det)) / (2 * ae);

      if (p2 > 0 && p2 < 1) {
        TA_tr = std::min(compute_AT_at_node(p2, indx, indy, indz), TA_tr);
      }
    }
  }

  return TA_tr;
}

SF_real eikonal_solver::solve_edges(mesh_int_t& indX, mesh_int_t& indY)
{
  SF::vector<SF_real> X(3), Y(3), XY(3);

  X[0] = x_coord[indX];
  X[1] = y_coord[indX];
  X[2] = z_coord[indX];
  Y[0] = x_coord[indY];
  Y[1] = y_coord[indY];
  Y[2] = z_coord[indY];

  XY[0] = Y[0] - X[0];
  XY[1] = Y[1] - X[1];
  XY[2] = Y[2] - X[2];

  Num  = mult_VtMV(XY, iD, XY);
  Num  = sqrt(Num);
  Frac = Num / CV_l;

  SF_real TAn = T_A[indY] + Frac;

  return TAn;
}  // end Solve_Edges

bool eikonal_solver::add_node_neighbor_to_list(SF_real& RT, SF_real& oldTA, SF_real& newTA)
{
  bool isNewTABetweenRTAndOldTA   = (RT < newTA) && (newTA < oldTA);
  bool isNewTABiggerRTBiggerOldTA = (oldTA < RT) && (RT < newTA);

  return (isNewTABetweenRTAndOldTA) || (isNewTABiggerRTBiggerOldTA);
}

SF_real eikonal_solver::mult_VtMV(SF::vector<SF_real>& X, SF::dmat<double>& iD, SF::vector<SF_real>& Y)
{
  SF_real result = 0;

  if (X[0] == Y[0] && X[1] == Y[1] && X[2] == Y[2]) {
    result = iD[0][0] * X[0] * X[0] + 2 * iD[0][1] * X[0] * X[1] + iD[1][1] * X[1] * X[1] + 2 * iD[0][2] * X[0] * X[2] + iD[2][2] * X[2] * X[2] + 2 * iD[1][2] * X[1] * X[2];
  } else {
    result = iD[0][0] * X[0] * Y[0] + iD[0][1] * (X[0] * Y[1] + X[1] * Y[0]) + iD[1][1] * X[1] * Y[1] + iD[0][2] * (X[0] * Y[2] + X[2] * Y[0]) + iD[2][2] * X[2] * Y[2] + iD[1][2] * (X[1] * X[2] + X[2] * Y[1]);
  }

  return result;

}  // mult_VtMV

SF_real eikonal_solver::compute_AT_at_node(SF_real& pp, mesh_int_t& indxf1, mesh_int_t& indyf1, mesh_int_t& indzf1)
{
  if (pp == 0) {
    Num  = sqrt(C);
    Frac = Num / CV_l;
    return T_A[indzf1] + Frac;

  } else if (pp == 1) {
    Num  = sqrt(A + B + C);
    Frac = Num / CV_l;
    return T_A[indyf1] + Frac;

  } else {
    Num    = sqrt(A * pp * pp + B * pp + C);
    Frac   = Num / CV_l;
    TA_ori = T_A[indyf1] * pp + T_A[indzf1] * (1 - pp);
    return TA_ori + Frac;
  }
}

SF_real eikonal_solver::compute_AT_at_node_in_tet(SF_real& pf2, SF_real& qf2, mesh_int_t indxf2, mesh_int_t indyf2, mesh_int_t indzf2, mesh_int_t indwf2)
{
  Num    = c1 * (pf2 * pf2) + c2 * (qf2 * qf2) + c3 * ((1 - pf2 - qf2) * (1 - pf2 - qf2)) + c4 * (pf2 * qf2) + c5 * (pf2 * (1 - pf2 - qf2)) + c6 * (qf2 * (1 - pf2 - qf2));
  Num    = sqrt(Num);
  Frac   = Num / CV_l;
  TA_ori = T_A[indyf2] * pf2 + T_A[indzf2] * qf2 + T_A[indwf2] * (1 - pf2 - qf2);
  return TA_ori + Frac;
}

void eikonal_solver::save_eikonal_state(const char* tsav_ext)
{
  if (get_rank() == 0) {
    FILE* file_writestate;
    char  buffer_writestate[1024];
    snprintf(buffer_writestate, sizeof buffer_writestate, "%s.%s.roe.dat", param_globals::write_statef, tsav_ext);
    file_writestate = fopen(buffer_writestate, "w");
    for (size_t jjj = 0; jjj < List.size(); jjj++) {
      fprintf(file_writestate, "%d %d %d %f %f %f %f \n", List[jjj], num_changes[jjj], nReadded2List[jjj], T_A[jjj], T_R[jjj], TA_old[jjj], D_I[jjj]);
    }

    fclose(file_writestate);
  }
}

void eikonal_solver::update_repolarization_times_from_rd(sf_vec& Vmv, sf_vec& Vmv_old, double time)
{
  SF_real* old_ptr = Vmv_old.ptr();
  SF_real* new_ptr = Vmv.ptr();

  for (int ind_nodes = 0; ind_nodes < num_pts; ind_nodes++) {
    double prev_R = T_R[ind_nodes];

    if (old_ptr[ind_nodes] >= param_globals::dream.repol_time_thresh && new_ptr[ind_nodes] < param_globals::dream.repol_time_thresh) {
      T_R[ind_nodes] = time;
    }
  }
  Vmv_old.release_ptr(old_ptr);
  Vmv.release_ptr(new_ptr);
}

void eikonal_solver::update_repolarization_times(const Ionics& ion)
{
  double t1, t0;
  get_time(t0);

  #ifdef _OPENMP
  int max_threads = omp_get_max_threads();
  omp_set_num_threads(1); // with the current setup of this function, it is better to run it in serial.
  #endif

  limpet::MULTI_IF* miif = ion.miif;

  const double time = user_globals::tm_manager->time,
               dt   = user_globals::tm_manager->time_step;

  for (int i = 0; i < miif->N_IIF; i++) {
    if (!miif->N_Nodes[i]) continue;
    limpet::IonIfBase* pIF     = ion.miif->IIF[i];
    limpet::IonIfBase* IIF_old = pIF->get_type().make_ion_if(pIF->get_target(),
                                                             pIF->get_num_node(),
                                                             ion.miif->plugtypes[i]);
    IIF_old->copy_SVs_from(*pIF, false);

    int current = 0;
    do {
      int ind_gb = (miif->NodeLists[i][current]);
      // Global index of current node;
      // save states at the moment
      double Vm_old   = miif->ldata[i][limpet::Vm][current];
      double Iion_old = miif->ldata[i][limpet::Iion][current];
      double prev_R   = T_R[ind_gb];

      bool cond_2upd = T_R[ind_gb] <= T_A[ind_gb] && T_A[ind_gb] != inf && Vm_old > param_globals::dream.repol_time_thresh;

      if (cond_2upd) {
        double elapsed_time = 0;
        double prev_Vm      = miif->ldata[i][limpet::Vm][current];

        do {
          elapsed_time += dt;
          pIF->compute(current, current + 1, miif->ldata[i]);
          miif->ldata[i][limpet::Vm][current] -= miif->ldata[i][limpet::Iion][current] * param_globals::dt;
          if (miif->ldata[i][limpet::Vm][current] < param_globals::dream.repol_time_thresh) {
            T_R[ind_gb] = time + elapsed_time;
            break;
          }
        } while (elapsed_time < 600);

        // Restore old states
        miif->ldata[i][limpet::Vm][current]   = Vm_old;
        miif->ldata[i][limpet::Iion][current] = Iion_old;
      }
      current++;
    } while (current < miif->N_Nodes[i]);
    pIF->copy_SVs_from(*IIF_old, false);
  }

  #ifdef _OPENMP
  omp_set_num_threads(max_threads); // restore max threads for parabolic solver
  #endif

  // treat solver statistics
  auto dur = timing(t1, t0);
  stats.slvtime_D += dur;
}

void eikonal_solver::compute_diffusion_current(const double& time, sf_vec& vm)
{
  double t1, t0;
  get_time(t0);

  SF_real* c = Idiff->ptr();
  SF_real* v = vm.ptr();
  SF_real  dT;
  double   e_on, e_off;

  int                           rank    = get_rank();
  const sf_mesh&                mesh    = get_mesh(intra_elec_msh);
  const SF::vector<mesh_int_t>& alg_nod = mesh.pl.algebraic_nodes();

  for (size_t j = 0; j < alg_nod.size(); j++) {
    mesh_int_t loc_nodal_idx = alg_nod[j];
    mesh_int_t loc_petsc_idx = local_nodal_to_local_petsc(mesh, rank, loc_nodal_idx);

    if ((time >= T_A[loc_nodal_idx]) && ((T_A[loc_nodal_idx] + 5) >= time) && (List[loc_nodal_idx] == 0)) {
      dT = time - T_A[loc_nodal_idx];

      switch (diff_cur[loc_nodal_idx].model) {
        case GAUSS:
          c[loc_petsc_idx] = diff_cur[loc_nodal_idx].alpha_1 * exp(-((dT - diff_cur[loc_nodal_idx].beta_1) / diff_cur[loc_nodal_idx].gamma_1) * ((dT - diff_cur[loc_nodal_idx].beta_1) / diff_cur[loc_nodal_idx].gamma_1))
                            + diff_cur[loc_nodal_idx].alpha_2 * exp(-((dT - diff_cur[loc_nodal_idx].beta_2) / diff_cur[loc_nodal_idx].gamma_2) * ((dT - diff_cur[loc_nodal_idx].beta_2) / diff_cur[loc_nodal_idx].gamma_2))
                            + diff_cur[loc_nodal_idx].alpha_3 * exp(-((dT - diff_cur[loc_nodal_idx].beta_3) / diff_cur[loc_nodal_idx].gamma_3) * ((dT - diff_cur[loc_nodal_idx].beta_3) / diff_cur[loc_nodal_idx].gamma_3));
          break;
        default:
          e_on             = (dT >= 0) ? 1.0 : 0.0;
          e_off            = (v[loc_petsc_idx] < diff_cur[loc_nodal_idx].V_th) ? 1.0 : 0.0;
          c[loc_petsc_idx] = diff_cur[loc_nodal_idx].A_F / diff_cur[loc_nodal_idx].tau_F * exp(dT / diff_cur[loc_nodal_idx].tau_F) * e_on * e_off;
      }
    }  // end if
  }  // end for loop

  Idiff->release_ptr(c);
  vm.release_ptr(v);

  // treat solver statistics
  auto dur = timing(t1, t0);
  stats.slvtime_B += dur;
}

void eikonal_solver::load_state_file()
{
  set_dir(INPUT);
  FILE* file_startstate;
  char  buffer_startstate[strlen(param_globals::start_statef) + 10];

  snprintf(buffer_startstate, sizeof buffer_startstate, "%s.dat", param_globals::start_statef);
  file_startstate = fopen(buffer_startstate, "r");

  if (file_startstate == NULL) {
    log_msg(NULL, 5, 0, "Not able to open state file: %s", buffer_startstate);
  } else if (param_globals::output_level) {
    log_msg(NULL, 0, 0, "Open state file for eikonal model: %s", buffer_startstate);
  }

  int   ListVal, numChangesVal, numChanges2Val;
  float T_AVal, T_RVal, TA_oldVal, D_IVal, PCLVal;
  int   index = 0;

  while (fscanf(file_startstate, "%d %d %d %f %f %f %f", &ListVal, &numChangesVal, &numChanges2Val, &T_AVal, &T_RVal, &TA_oldVal, &D_IVal) == 7) {
    List[index]          = ListVal;
    num_changes[index]   = numChangesVal;
    nReadded2List[index] = numChanges2Val;
    T_A[index]           = T_AVal;
    T_R[index]           = T_RVal;
    TA_old[index]        = TA_oldVal;
    D_I[index]           = D_IVal;
    ++index;
  }
  if (param_globals::output_level) log_msg(NULL, 0, 0, "Number of nodes in active list: %i", sum(List));

  fclose(file_startstate);
}

void eikonal_solver::init_diffusion_current()
{
  const sf_mesh&         mesh = get_mesh(intra_elec_msh);
  SF::vector<mesh_int_t> vertices;

  for (int num_i = 0; num_i < param_globals::num_imp_regions; num_i++) {
    for (int num_id = 0; num_id <= param_globals::imp_region[num_i].num_IDs; num_id++) {
      indices_from_region_tag(vertices, mesh, param_globals::imp_region[num_i].ID[num_id]);
      for (int i = 0; i < vertices.size(); i++) {
        diff_cur[vertices[i]].model = static_cast<eikonal_solver::Idiff_t>(param_globals::imp_region[num_i].dream.Idiff.model);
        switch (diff_cur[vertices[i]].model) {
          case GAUSS:
            diff_cur[vertices[i]].alpha_1 = param_globals::imp_region[num_i].dream.Idiff.alpha_i[0];
            diff_cur[vertices[i]].alpha_2 = param_globals::imp_region[num_i].dream.Idiff.alpha_i[1];
            diff_cur[vertices[i]].alpha_3 = param_globals::imp_region[num_i].dream.Idiff.alpha_i[2];
            diff_cur[vertices[i]].beta_1  = param_globals::imp_region[num_i].dream.Idiff.beta_i[0];
            diff_cur[vertices[i]].beta_2  = param_globals::imp_region[num_i].dream.Idiff.beta_i[1];
            diff_cur[vertices[i]].beta_3  = param_globals::imp_region[num_i].dream.Idiff.beta_i[2];
            diff_cur[vertices[i]].gamma_1 = param_globals::imp_region[num_i].dream.Idiff.gamma_i[0];
            diff_cur[vertices[i]].gamma_2 = param_globals::imp_region[num_i].dream.Idiff.gamma_i[1];
            diff_cur[vertices[i]].gamma_3 = param_globals::imp_region[num_i].dream.Idiff.gamma_i[2];
            break;
          default:
            diff_cur[vertices[i]].A_F   = param_globals::imp_region[num_i].dream.Idiff.A_F;
            diff_cur[vertices[i]].tau_F = param_globals::imp_region[num_i].dream.Idiff.tau_F;
            diff_cur[vertices[i]].V_th  = param_globals::imp_region[num_i].dream.Idiff.V_th;
        };
      }
    }
  }
}

void eikonal_solver::translate_stim_to_eikonal()
{
  // Set stimulus for eikonal model
  if (param_globals::output_level) log_msg(NULL, 0, 0, "Initializing stimuli for eikonal model ...");

  SF::vector<SF_real>    StartStimulus;
  SF::vector<mesh_int_t> IndexStimulus;
  int                    number_stimulus = param_globals::num_stim;
  int                    total_stimuli   = 0;
  int                    count_stm       = 0;
  int                    count           = 0;
  Index_currStim                         = 0;

  for (stimulus s : *stimuliRef) {
    total_stimuli += s.ptcl.npls;
    for (int idx_pls = 0; idx_pls < s.ptcl.npls; idx_pls++) {
      SF_real start_time = s.ptcl.start + s.ptcl.pcl * idx_pls;
      StartStimulus.push_back(start_time);
      IndexStimulus.push_back(s.idx);
      count_stm++;
    }
  }
  // sort by start times in ascending order
  bubble_sort(StartStimulus, IndexStimulus, total_stimuli);

  for (int ind_stim = 0; ind_stim < total_stimuli; ind_stim++) {
    stimulus s = (*stimuliRef)[IndexStimulus[ind_stim]];
    // Loop over vertices for each defined stim
    for (int ind_nodes = 0; ind_nodes < s.electrode.vertices.size(); ind_nodes++) {
      StimulusPoints[count] = s.electrode.vertices[ind_nodes];
      StimulusTimes[count]  = StartStimulus[ind_stim];
      count++;
    }
  }
  // Add -1 at the end to indicate we ran out of stimulations
  StimulusPoints.resize(count + 2, -1);  // Nodes that have initial stimulus
  StimulusTimes.resize(count + 2, -1);   // Times when the initial nodes are stimulated
}

void eikonal_solver::create_node_to_node_graph()
{
  // Create the Nodes to Nodes Graph (Per Each node the graph finds the NB nodes)
  n2n_dsp.resize(num_pts + 1, 0);
  const sf_mesh& mesh = get_mesh(intra_elec_msh);
  twoFib = mesh.she.size() > 0;

  for (int point_idx = 0; point_idx < num_pts; point_idx++) {
    int                     numNBElem = n2e_dsp[point_idx + 1] - n2e_dsp[point_idx];
    std::vector<mesh_int_t> NNNodesxNodes(numNBElem * numPtsxElem);

    for (int eedsp = 0; eedsp < numNBElem; eedsp++) {
      int currElem = n2e_con[n2e_dsp[point_idx] + eedsp];

      if (twoFib) {
        fiber_node[3 * point_idx]     = mesh.fib[3 * currElem];
        fiber_node[3 * point_idx + 1] = mesh.fib[3 * currElem + 1];
        fiber_node[3 * point_idx + 2] = mesh.fib[3 * currElem + 2];

        sheet_node[3 * point_idx]     = mesh.she[3 * currElem];
        sheet_node[3 * point_idx + 1] = mesh.she[3 * currElem + 1];
        sheet_node[3 * point_idx + 2] = mesh.she[3 * currElem + 2];

      } else {
        fiber_node[3 * point_idx]     = mesh.fib[3 * currElem];
        fiber_node[3 * point_idx + 1] = mesh.fib[3 * currElem + 1];
        fiber_node[3 * point_idx + 2] = mesh.fib[3 * currElem + 2];
      }

      for (int nndsp = 0; nndsp < numPtsxElem; nndsp++) {
        NNNodesxNodes[eedsp * numPtsxElem + nndsp] = e2n_con[elem_start[currElem] + nndsp];
      }
    }

    std::sort(NNNodesxNodes.begin(), NNNodesxNodes.end());
    auto it = std::unique(NNNodesxNodes.begin(), NNNodesxNodes.end());
    NNNodesxNodes.erase(it, NNNodesxNodes.end());
    it = std::remove(NNNodesxNodes.begin(), NNNodesxNodes.end(), point_idx);
    NNNodesxNodes.erase(it, NNNodesxNodes.end());
    int max = NNNodesxNodes.size();

    n2n_dsp[point_idx + 1] = n2n_dsp[point_idx] + NNNodesxNodes.size();

    n2n_connect.insert(n2n_connect.end(), NNNodesxNodes.begin(), NNNodesxNodes.end());
  }

  n2n_dsp[num_pts] = n2n_connect.size();
}

void node_stats::init_logger(const char* filename)
{
  logger = f_open(filename, "w");

  const char* h1 = "          ------ ---------- ---------- ------- ------- | List logic ----- -------- | Neighbor node ---- |";
  const char* h2 = "           cycle         AT     old AT      RT      DI |  Status    Entry     Exit |      ID         AT |";

  if (logger == NULL)
    log_msg(NULL, 3, 0, "%s error: Could not open file %s in %s. Turning off logging.\n",
            __func__, filename);
  else {
    log_msg(logger, 0, 0, "%s", h1);
    log_msg(logger, 0, 0, "%s", h2);
  }
}

void node_stats::log_stats(double time, bool cflg)
{
  if (!this->logger) return;

  char abuf[256];
  char bbuf[256];
  char cbuf[256];

  // create nicer output in logger for inf, -inf, nan
  std::ostringstream oss_TA, oss_TA_, oss_TR, oss_DI, oss_nbnTA;
  oss_TA << this->T_A;
  oss_TA_ << this->T_A_;
  oss_TR << this->T_R;
  oss_DI << this->D_I;
  oss_nbnTA << this->nbn_T_A;

  if (this->idXNB == std::numeric_limits<Int>::min()) {
    snprintf(cbuf, sizeof cbuf, "%7s %10s", "-", "-");
  } else {
    snprintf(cbuf, sizeof cbuf, "%7d %10s", this->idXNB, oss_nbnTA.str().c_str());
  }

  snprintf(abuf, sizeof abuf, "%6d %10s %10s %7s %7s", this->cycle, oss_TA.str().c_str(), oss_TA_.str().c_str(), oss_TR.str().c_str(), oss_DI.str().c_str());
  snprintf(bbuf, sizeof bbuf, "%7s %8s %8s", this->status, this->reasonIn, this->reasonOut);

  unsigned char flag = cflg ? ECHO : 0;
  log_msg(this->logger, 0, flag | FLUSH | NONL, "%9.3f %s | %s | %s |\n", time, abuf, bbuf, cbuf);

  this->reasonIn  = "-";
  this->reasonOut = "-";
  this->T_A_      = this->T_A;
  this->T_A       = std::numeric_limits<double>::quiet_NaN();
  this->T_R       = std::numeric_limits<double>::quiet_NaN();
  this->D_I       = std::numeric_limits<double>::quiet_NaN();
  this->idXNB     = std::numeric_limits<Int>::min();
  this->nbn_T_A   = std::numeric_limits<double>::quiet_NaN();
  this->cycle++;
}

void node_stats::update_status(enum status s, enum reason r)
{
  // directly convert enums to strings for logging
  const char* stat_str;
  const char* reas_str;
  switch (s) {
    case node_stats::in: stat_str = "in"; break;
    case node_stats::out: stat_str = "out"; break;
  }

  switch (r) {
    case node_stats::none: reas_str = "-"; break;
    case node_stats::nbn: reas_str = "nbn"; break;
    case node_stats::conv: reas_str = "conv"; break;
    case node_stats::stim: reas_str = "stim"; break;
  }

  this->status = stat_str;
  if (s == in) {
    this->reasonIn = reas_str;
  } else {
    this->reasonOut = reas_str;
  }
}

}  // namespace opencarp