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 "electric_integrators.h"
#include <cstdlib>
#include "SF_globals.h"
#include "basics.h"
#include "petsc_utils.h"
#include "timers.h"
#include "stimulate.h"
 
#include "SF_init.h" 
 
namespace opencarp
{
 
void Eikonal::initialize()
{
  if (get_size() > 1) {
    log_msg(NULL, 5, 0, "DREAM/Eikonal physics currently do not support MPI parallelization. Use openMP instead. Aborting!");
    EXIT(EXIT_FAILURE);
  }
  
  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. The standalone eikonal solver does not require 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();
 
  this->initialize_time += timing(t2, t1);
  log_msg(NULL, 0, 0, "All done in %f sec.", float(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;
  }
 
  {
    mesh_t      mt   = g == Eikonal::intra_grid ? intra_elec_msh : extra_elec_msh;
    const char* file = g == Eikonal::intra_grid ? param_globals::gi_scale_vec : param_globals::ge_scale_vec;
    if (strlen(file))
      read_el_scale_vec(file, mt, m->el_scale, m->el_scale_dpn);
  }
}
 
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()
{
  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();
      // destroy ionics before closing the logger: the ionic interface holds an alias of it
      ion.destroy();
  }
 
  // close logger
  f_close(logger);
}
 
void Eikonal::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 (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);
 
    if (s.electrode.dump_vtx) 
      s.dump_vtx_file(i);
 
    log_msg(NULL, 2, 0, "Only geometry, start time, npls, and bcl of stim[%i] are used", i);
 
    if (param_globals::stim[i].pulse.dumpTrace && get_rank() == 0) {
      set_dir(OUTPUT);
      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();
      break;
    default:
      parab_solver.init();
      parab_solver.rebuild_matrices(mtype, *ion.miif, logger);
      eik_solver.init();
      eik_solver.init_imp_region_properties();
      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::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");
  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_%lld.dat", static_cast<long long>(nodeData.idX));
    nodeData.init_logger(buf);
  }
 
  // currently only used for output
  const sf_mesh& mesh = get_mesh(intra_elec_msh);
  twoFib              = mesh.she.size() > 0;
  SF::init_vector(&Idiff, mesh, 1, sf_vec::algebraic);
  SF::init_vector(&AT, mesh, 1, sf_vec::algebraic);
  SF::init_vector(&RT, mesh, 1, 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);
 
  List.assign(num_pts, 0);           // List of active nodes
  T_A.assign(num_pts, inf);          // Activation times
  D_I.assign(num_pts, 0);            // Diastolic Intervals
  T_R.assign(num_pts, -400);         // Recovery times
  TA_old.assign(num_pts, inf);       // Activation Time is the previous cycle
  num_changes.assign(num_pts, 0);    // Number of Times entered in the List per current LAT
  nReadded2List.assign(num_pts, 0);  // Number of Times entered in the List per current LAT
  stim_status.assign(num_pts, 0);    // Status of nodes
 
  switch (mesh.type[0]) {
    case 0:
      MESH_SIZE = 4;
      break;
 
    case 6:
      MESH_SIZE = 3;
      break;
 
    case 7:
      MESH_SIZE = 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();
 
  precompute_squared_anisotropy_metric();
}
 
void eikonal_solver::init_imp_region_properties()
{
  double t1, t2;
  get_time(t1);
 
  const sf_mesh& mesh = get_mesh(intra_elec_msh);
 
  diff_cur.resize(mesh.l_numpts);
  rho_cvrest.resize(mesh.l_numpts);
  kappa_cvrest.resize(mesh.l_numpts);
  theta_cvrest.resize(mesh.l_numpts);
  denom_cvrest.resize(mesh.l_numpts);
 
  // --- Precompute resolved region IDs for each point ---
  // Use the same region delegation as the ionics class
  SF::vector<int>         regionIDs;
  SF::vector<RegionSpecs> rs(param_globals::num_imp_regions);
  for (size_t i = 0; i < rs.size(); i++) {
    rs[i].nsubregs   = param_globals::imp_region[i].num_IDs;
    rs[i].subregtags = param_globals::imp_region[i].ID;
    for (int j = 0; j < rs[i].nsubregs; j++) {
      if (rs[i].subregtags[j] == -1 && get_rank() == 0)
        log_msg(NULL, 3, ECHO, "Warning: not all %u IDs provided for imp_region[%u]!\n", rs[i].nsubregs, i);
    }
  }
  if (rs.size() == 1) {
    regionIDs.assign(mesh.l_numpts, 0);
  } else {
    region_mask(intra_elec_msh, rs, regionIDs, false, "imp_regions");
  }
  
  // --- Parallel initialization, one write per vertex ---
  #pragma omp parallel for schedule(dynamic)
  for (int v = 0; v < mesh.l_numpts; v++) {
    int reg = regionIDs[v];
 
    const auto& region_diff = param_globals::imp_region[reg].dream.Idiff;
    const auto& region_rest = param_globals::imp_region[reg].dream.CVrest;
 
    auto model        = static_cast<eikonal_solver::Idiff_t>(region_diff.model);
    diff_cur[v].model = model;
 
    if (model == GAUSS) {
      diff_cur[v].alpha_1 = region_diff.alpha_i[0];
      diff_cur[v].alpha_2 = region_diff.alpha_i[1];
      diff_cur[v].alpha_3 = region_diff.alpha_i[2];
      diff_cur[v].beta_1  = region_diff.beta_i[0];
      diff_cur[v].beta_2  = region_diff.beta_i[1];
      diff_cur[v].beta_3  = region_diff.beta_i[2];
      diff_cur[v].gamma_1 = region_diff.gamma_i[0];
      diff_cur[v].gamma_2 = region_diff.gamma_i[1];
      diff_cur[v].gamma_3 = region_diff.gamma_i[2];
    } else {
      diff_cur[v].A_F   = region_diff.A_F;
      diff_cur[v].tau_F = region_diff.tau_F;
      diff_cur[v].V_th  = region_diff.V_th;
    }
 
    rho_cvrest[v]   = region_rest.rho;
    kappa_cvrest[v] = region_rest.kappa;
    theta_cvrest[v] = region_rest.theta;
    denom_cvrest[v] = log(region_rest.rho) / region_rest.psi;
  }
  if (param_globals::output_level) log_msg(NULL, 0, 0, "Diffusion current and CV restitution initialized in %f sec.", timing(t2, t1));
}
 
void eikonal_solver::translate_stim_to_eikonal()
{
    double t1, t2;
    get_time(t1);
    Index_currStim = 0;
 
    // Collect all pulses as (start_time, stimulus_index)
    std::vector<std::pair<SF_real, int>> stim_events;
    stim_events.reserve(param_globals::num_stim * 8); // rough guess
 
    for (int stim_idx = 0; stim_idx < stimuliRef->size(); ++stim_idx) {
        const stimulus& s = (*stimuliRef)[stim_idx];
        for (int idx_pls = 0; idx_pls < s.ptcl.npls; ++idx_pls) {
            SF_real start_time = s.ptcl.start + s.ptcl.pcl * idx_pls;
            stim_events.emplace_back(start_time, stim_idx);
        }
    }
 
    // Sort events by start_time
    std::sort(stim_events.begin(), stim_events.end(),
              [](auto& a, auto& b) { return a.first < b.first; });
 
    // Reserve enough space for output
    size_t total_nodes = 0;
    for (auto& ev : stim_events)
        total_nodes += (*stimuliRef)[ev.second].electrode.vertices.size();
    StimulusPoints.resize(total_nodes + 1, -1);
    StimulusTimes.resize(total_nodes + 1, -1);
 
    // Fill outputs
    size_t count = 0;
    for (auto& ev : stim_events) {
        const stimulus& s = (*stimuliRef)[ev.second];
        for (mesh_int_t v : s.electrode.vertices) {
            StimulusPoints[count] = v;
            StimulusTimes[count]  = ev.first;
            ++count;
        }
    }
 
    if (param_globals::output_level)
        log_msg(NULL, 0, 0, "Translating stimuli for eikonal model done in %f sec.", timing(t2, t1));
}
 
void eikonal_solver::create_node_to_node_graph()
{
  double t1, t2;
  get_time(t1);
 
  n2n_dsp.resize(num_pts + 1, 0);
  const sf_mesh& mesh = get_mesh(intra_elec_msh);
 
  // Thread-local storage for neighbors
  std::vector<std::vector<mesh_int_t>> thread_neighbors(num_pts);
 
  // Preallocate thread-local marker arrays for O(n) duplicate removal
  #pragma omp parallel
  {
    std::vector<char> mark(num_pts, 0);  // one per thread
 
    #pragma omp for schedule(dynamic)
    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> neighbors;
      neighbors.reserve(numNBElem * MESH_SIZE);  // rough estimate
 
      // Loop over all elements containing this node
      for (int eedsp = 0; eedsp < numNBElem; eedsp++) {
        int currElem = n2e_con[n2e_dsp[point_idx] + eedsp];
 
        // Loop over all nodes of the current element
        for (int nndsp = 0; nndsp < MESH_SIZE; nndsp++) {
          int nb = e2n_con[elem_start[currElem] + nndsp];
 
          if (nb == point_idx) continue;  // skip self
          if (!mark[nb]) {
            mark[nb] = 1;
            neighbors.push_back(nb);
          }
        }
      }
 
      // Reset marks for the next iteration
      for (int nb : neighbors) mark[nb] = 0;
 
      // Move neighbor list into the final container
      thread_neighbors[point_idx] = std::move(neighbors);
    }
  }
 
  // Build prefix sums (serial, cheap compared to above)
  for (int i = 0; i < num_pts; i++) {
    n2n_dsp[i + 1] = n2n_dsp[i] + thread_neighbors[i].size();
  }
 
  // Allocate once, then fill in parallel
  n2n_connect.resize(n2n_dsp[num_pts]);
 
  #pragma omp parallel for schedule(static)
  for (int i = 0; i < num_pts; i++) {
    std::copy(thread_neighbors[i].begin(),
              thread_neighbors[i].end(),
              n2n_connect.begin() + n2n_dsp[i]);
  }
 
  if (param_globals::output_level) {
    log_msg(NULL, 0, 0, "Node-to-node graph done in %f sec.", timing(t2, t1));
  }
}
 
void eikonal_solver::precompute_squared_anisotropy_metric()
{
  double t1, t2;
  get_time(t1);
 
  const sf_mesh& mesh = get_mesh(intra_elec_msh);
 
  S.resize(mesh.l_numelem);     // store anisotropy matrices (dimensionless)
  CV_L.resize(mesh.l_numelem);  // store CV_L per element for later use
 
  // temporary variables
  SF::dmat<double> I(3, 3);
  I.assign(0.0);
  I.diag(1.0);
 
  #pragma omp parallel
  {
    SF::dmat<double> Saniso(3, 3);
    SF::dmat<double> Ff(3, 3);
    SF::dmat<double> diff(3, 3);
    SF::Point        f, s, n;
 
    #pragma omp for schedule(static)
    for (int eidx = 0; eidx < mesh.l_numelem; eidx++) {
      double vl = 0.0, vt = 0.0, vn = 0.0;
 
      // Find region match (break early when found)
      for (int g = 0; g < param_globals::num_gregions; g++) {
        const auto& reg   = param_globals::gregion[g];
        const auto& dream = reg.dream;
        for (int j = 0; j < reg.num_IDs; j++) {
          if (mesh.tag[eidx] == reg.ID[j]) {
            vl = dream.vel_l;
            vt = dream.vel_t;
            vn = dream.vel_n;
            goto region_found;  // break out of both loops
          }
        }
      }
      region_found:;
 
      // Store CV_L directly (for later rescaling)
      CV_L[eidx] = vl;
 
      // Compute anisotropy ratios relative to CV_L
      double AR_T2 = (vl / vt) * (vl / vt);
      double AR_N2 = (vl / vn) * (vl / vn);
 
      f.x = mesh.fib[3 * eidx + 0];
      f.y = mesh.fib[3 * eidx + 1];
      f.z = mesh.fib[3 * eidx + 2];
 
      if (twoFib) {
        s.x = mesh.she[3 * eidx + 0];
        s.y = mesh.she[3 * eidx + 1];
        s.z = mesh.she[3 * eidx + 2];
        n   = cross(f, s);
 
        SF::outer_prod(f, f, 1.0, Saniso[0], false);
        SF::outer_prod(s, s, AR_T2, Saniso[0], true);
        SF::outer_prod(n, n, AR_N2, Saniso[0], true);
      } else {
        SF::outer_prod(f, f, 1.0, Ff[0], false);
        diff = I - Ff;
        diff *= AR_T2;
 
        Saniso = Ff + diff;
      }
 
      S[eidx] = Saniso;
    }
  }
 
  if (param_globals::output_level)
    log_msg(NULL, 0, 0, "Anisotropy tensors precomputed in %f sec.", timing(t2, t1));
}
 
void eikonal_solver::FIM()
{
  double t0, t1;
  get_time(t0);
 
  // --- 0) Reset all activation times to infinity
  std::fill(T_A.begin(), T_A.end(), inf);
 
  // --- 1) Build initial active list = neighbors of stimuli
  std::vector<mesh_int_t> activeList;
  activeList.reserve(num_pts / 10);  // heuristic reserve to avoid frequent reallocs
  std::vector<char> in_active(num_pts, 0);
 
  // Mark stimuli 
  for (size_t si = 0; si < StimulusTimes.size() - 1; ++si) { // skip last entry, since its a -1 placeholder currently still used in DREAM
    mesh_int_t s = StimulusPoints[si];
    T_A[s]       = StimulusTimes[si];
  }
  // Seed neighbors; separate loop to avoid adding stim nodes into the list if they are a neighbor
  for (size_t si = 0; si < StimulusTimes.size() - 1; ++si) { // skip last entry, since its a -1 placeholder currently still used in DREAM
    mesh_int_t s = StimulusPoints[si];
    for (int off = n2n_dsp[s]; off < n2n_dsp[s + 1]; ++off) {
      mesh_int_t nb = n2n_connect[off];
      if (T_A[nb] == inf) {
        add_to_active(activeList, in_active, nb);
      }
    }
  }
 
  // --- 2) Iterative solve of activeList
  int niter = 0;
  while (!activeList.empty() && niter <= param_globals::dream.fim.max_iter) {
    const std::vector<mesh_int_t> activeVec = std::move(activeList);
    activeList.clear();
 
    #pragma omp parallel
    {
      std::vector<mesh_int_t> local_active;
      local_active.reserve(64);  // small local buffer
 
      #pragma omp for schedule(dynamic)
      for (size_t i = 0; i < activeVec.size(); ++i) {
        mesh_int_t id = activeVec[i];
        remove_from_active(in_active, id);
 
        SF_real p = T_A[id];
        SF_real q = update(id);
 
        #pragma omp atomic write
        T_A[id] = q;
 
        if (std::fabs(p - q) < param_globals::dream.fim.tol) {
          for (int off = n2n_dsp[id]; off < n2n_dsp[id + 1]; ++off) {
            mesh_int_t nb = n2n_connect[off];
            p             = T_A[nb];
            q             = update(nb);
            if (p > q) {
              #pragma omp atomic write
              T_A[nb] = q;
              local_active.push_back(nb);
            }
          }
        } else {
          local_active.push_back(id);  // non-converged -> add back
        }
      }
 
      // Merge local results
      #pragma omp critical
      {
        for (mesh_int_t nb : local_active) {
          add_to_active(activeList, in_active, nb);
        }
      }
    }
    ++niter;
  }
 
  // collect iteration stats
  stats.update_iter(niter);
  auto [minIt, maxIt] = std::minmax_element(T_A.begin(), T_A.end());
  actMIN              = *minIt;
  actMAX              = *maxIt;
 
  // --- 3) Copy into AT for output
  double*        atc = AT->ptr();
  const SF_real* t_a = T_A.data();
  #pragma omp parallel for simd
  for (mesh_int_t i = 0; i < num_pts; ++i)
    atc[i] = (t_a[i] == inf ? -1.0 : t_a[i]);
  AT->release_ptr(atc);
 
  // --- 4) Timing & list‐size stats
  auto dur = timing(t1, t0);
  stats.slvtime_A += dur;
  stats.minAT      = actMIN;
  stats.maxAT      = actMAX;
  stats.activeList = activeList.size();
  stats.bc_status  = true;
  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();
  }
 
  double 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::update(mesh_int_t& indX, SF_real CVrest_factor, bool isDREAM)
{
  switch (MESH_SIZE) {
    case 4: return update_impl<4>(indX, CVrest_factor, isDREAM);
    case 3: return update_impl<3>(indX, CVrest_factor, isDREAM);
    case 2: return update_impl<2>(indX, CVrest_factor, isDREAM);
    default:
      return T_A[indX]; // fallback
  }
}
 
template <int N>
SF_real eikonal_solver::update_impl(mesh_int_t& indX, SF_real CVrest_factor, bool CheckValidity)
{
  double         min  = inf;
  const double   time = user_globals::tm_manager->time;
  const sf_mesh& mesh = get_mesh(intra_elec_msh);
 
  // element-wise update
  for (int e_i = n2e_dsp[indX]; e_i < n2e_dsp[indX + 1]; e_i++) {  // gather all elements the node belongs to
    int                                        Elem_i = n2e_con[e_i];
    mesh_int_t                                 indEle = elem_start[Elem_i];
    std::array<SF::Point, N>                   base;    // points
    std::array<double, N>                      values;  // activation times
    std::array<int, N>                         nodeIDs;
 
    std::size_t                                k = 0;
    for (std::size_t j = 0; j < N; j++) {  // loop over nodes of element e_i
      int n_i = e2n_con[indEle + j];
      if (n_i != indX) {
        base[k].x = mesh.xyz[3 * n_i + 0]; base[k].y = mesh.xyz[3 * n_i + 1]; base[k].z = mesh.xyz[3 * n_i + 2];
        values[k] = T_A[n_i];
        nodeIDs[k] = n_i;
        k++;
      } else {  // last slot is reserved for the current vertex we are solving for
        base[N - 1].x = mesh.xyz[3 * indX + 0]; base[N - 1].y = mesh.xyz[3 * indX + 1]; base[N - 1].z = mesh.xyz[3 * indX + 2];
        values[N - 1] = T_A[indX];
        nodeIDs[N - 1] = indX;
      }
    }
    // scale slowness metric by CV
    double cv = CV_L[Elem_i] * CVrest_factor;
    if (cv == 0.0) continue; // avoid 1/(cv*cv)
    SF::dmat<double> D = 1.0 / (cv * cv) * S[Elem_i];
 
    // 1) solve full N-simplex (eikonal run OR if valid for DREAM)
    if (!CheckValidity || !is_not_valid_update<N>(nodeIDs, time)) {
      LocalSolver<N> solver(D, base, values);
      SF_real tmp = solver.solve();
      if (min > tmp && tmp > T_R[indX] && compute_H(indX, tmp) > 0.0) {
        min = tmp;
      }
      continue;
    }
 
    // If the full N-simplex was not valid, we have to do additional checks for the DREAM,
    // since a subsimplex could be valid if e.g. only one node of a tet/triangle is invalid
    // 2) triangles that include the target (only done if MESH_SIZE is 4)
    if constexpr (N - 1 == 3) {
      // neighbor indices are 0..(N-2); choose pairs (i,j) and add target (N-1)
      for (int i = 0; i < (N - 1); ++i) {
        for (int j = i + 1; j < (N - 1); ++j) {
          // build nodeIDs/points/values for face {i, j, target}
          std::array<int, 3> tri_ids{nodeIDs[i], nodeIDs[j], nodeIDs[N - 1]};
          if (is_not_valid_update<3>(tri_ids, time)) continue;
 
          std::array<SF::Point, 3> tri_pts{base[i], base[j], base[N - 1]};
          std::array<double, 3>    tri_vals{values[i], values[j], values[N - 1]};
 
          LocalSolver<3> solver(D, tri_pts, tri_vals);
          // min = std::min(min, solver.solve());
          SF_real tmp = solver.solve();
          if (min > tmp && tmp > T_R[indX] && compute_H(indX, tmp) > 0.0) {
            min = tmp;
          }
        }
      }
    }
 
    // 3) edges that include the target
    for (int i = 0; i < N - 1; i++) {
      std::array<int, 2> edge_ids{nodeIDs[i], nodeIDs[N - 1]};
      if (is_not_valid_update<2>(edge_ids, time)) continue;
 
      std::array<SF::Point, 2> edge_pts{base[i], base[N - 1]};
      std::array<double, 2>    edge_vals{values[i], values[N - 1]};
 
      LocalSolver<2> solver(D, edge_pts, edge_vals);
      SF_real        tmp = solver.solve();
      if (min > tmp && tmp > T_R[indX] && compute_H(indX, tmp) > 0.0) {
        min = tmp;
      }
    }
  }
 
  return min;
}
 
template <int N>
bool eikonal_solver::is_not_valid_update(const std::array<int, N>& nodeIDs, double time)
{
  const int  target     = nodeIDs[N - 1];
  const bool failedStim = (stim_status[target] == 2);
 
  for (int i = 0; i < N - 1; i++) {
    const int nb = nodeIDs[i];
 
    if (T_A[nb] < time) return true;
    if (T_A[nb] < T_R[nb]) return true;
    if (failedStim && stim_status[nb] == 1) return true;
  }
 
  return false;
}
 
SF_real eikonal_solver::compute_H(mesh_int_t& indX, SF_real& tmpTA)
{
  // Apply a refractory delay if stimulus previously failed
  const double delay = (stim_status[indX] == 2) ? 5.0 : 0.0;
 
  // Early exit conditions
  if (tmpTA == inf ||
      (T_R[indX] + delay) == -400.0 ||
      std::fabs(tmpTA - TA_old[indX]) < param_globals::dream.fim.tol ||
      T_R[indX] > tmpTA) {
    return 1.0;
  }
 
  // Compute diastolic interval
  const double DI = tmpTA - (T_R[indX] + delay);
  D_I[indX]       = DI;
 
  // CV restitution factor
  const double exponent  = -(DI + kappa_cvrest[indX]) * denom_cvrest[indX];
  const double Factor_DI = 1.0 - rho_cvrest[indX] * std::exp(exponent);
 
  if (Factor_DI <= 0.0) {
    return 0.0;
  }
 
  // Apply restitution curve behavior
  return (DI <= theta_cvrest[indX]) ? -Factor_DI : Factor_DI;
}
 
SF_real eikonal_solver::compute_coherence(mesh_int_t& indX)
{
  SF_real scaling = 1.0;
  SF_real p       = T_A[indX];                    // starting guess (previous arrival time)
  SF_real q       = update(indX, scaling, true);  // candidate new arrival time
 
  for (int niter = 0; niter < param_globals::dream.fim.max_coh; ++niter) {
    scaling = compute_H(indX, p);                      // restitution-based scaling for CV (can be negative mid-iteration)
    q       = update(indX, std::fabs(scaling), true);  // new candidate arrival time; use fabs() to allow intermediate negative scaling
 
    if (std::fabs(p - q) < param_globals::dream.fim.tol)
      break;  // converged
 
    p = q;  // continue iterating
  }
 
  // Only accept q if final state is physiologically valid
  return (compute_H(indX, q) < 0.0) ? T_A[indX] : 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) {
        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
 
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::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;
    }
  }
}
 
bool eikonal_solver::add_node_neighbor_to_list(SF_real& RT, SF_real& oldTA, SF_real& newTA)
{
  // Condition (A): RT < newTA < oldTA
  // Condition (B): oldTA < RT < newTA
  return ((RT < newTA) && (newTA < oldTA)) ||
         ((oldTA < RT) && (RT < newTA));
}
 
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();
 
  const double thresh = param_globals::dream.repol_time_thresh;
 
  for (int ind_nodes = 0; ind_nodes < num_pts; ind_nodes++) {
    if (old_ptr[ind_nodes] >= thresh && new_ptr[ind_nodes] < 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 t0, t1;
  get_time(t0);
 
  SF_real* c = Idiff->ptr();
  SF_real* v = vm.ptr();
 
  const sf_mesh&                mesh    = get_mesh(intra_elec_msh);
  const SF::vector<mesh_int_t>& alg_nod = mesh.pl.algebraic_nodes();
  int                           rank    = get_rank();
 
  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);
 
    double TA = T_A[loc_nodal_idx];
    if (!(time >= TA && (TA + 5) >= time && List[loc_nodal_idx] == 0))
      continue;
 
    double dT   = time - TA;
    auto&  node = diff_cur[loc_nodal_idx];
 
    switch (node.model) {
      case GAUSS: {
        double term1     = (dT - node.beta_1) / node.gamma_1;
        double term2     = (dT - node.beta_2) / node.gamma_2;
        double term3     = (dT - node.beta_3) / node.gamma_3;
        c[loc_petsc_idx] = node.alpha_1 * exp(-term1 * term1) + node.alpha_2 * exp(-term2 * term2) + node.alpha_3 * exp(-term3 * term3);
        break;
      }
      default: {
        double e_on      = (dT >= 0.0) ? 1.0 : 0.0;
        double e_off     = (v[loc_petsc_idx] < node.V_th) ? 1.0 : 0.0;
        c[loc_petsc_idx] = node.A_F / node.tau_F * exp(dT / node.tau_F) * e_on * e_off;
        break;
      }
    }
  }
 
  Idiff->release_ptr(c);
  vm.release_ptr(v);
 
  stats.slvtime_B += timing(t1, t0);
}
 
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, "%lld %lld %lld %f %f %f %f \n",
              static_cast<long long>(List[jjj]),
              static_cast<long long>(num_changes[jjj]),
              static_cast<long long>(nReadded2List[jjj]),
              T_A[jjj], T_R[jjj], TA_old[jjj], D_I[jjj]);
    }
 
    fclose(file_writestate);
  }
}
 
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 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, "%7lld %10s", static_cast<long long>(this->idXNB), oss_nbnTA.str().c_str());
  }
 
  snprintf(abuf, sizeof abuf, "%6lld %10s %10s %7s %7s",
           static_cast<long long>(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;
  }
}
 
template<int MESH_SIZE>
double LocalSolver<MESH_SIZE>::solve() {
  const SF::Point& x1 = points[0];
  const SF::Point& x2 = points[1];
 
  const double& u1 = values[0];
 
  if constexpr (MESH_SIZE == 2) {
    return tsitsiklis_update_line({x1,x2}, D, u1);
 
  } else if constexpr (MESH_SIZE == 3) {
    const SF::Point& x3 = points[2];
    const double&    u2 = values[1];
    return tsitsiklis_update_triangle({x1,x2,x3}, D, {u1,u2});
 
  } else if constexpr (MESH_SIZE == 4) {
    const SF::Point& x3 = points[2];
    const SF::Point& x4 = points[3];
    const double&    u2 = values[1];
    const double&    u3 = values[2];
 
    double u_tet = tsitsiklis_update_tetra({x1,x2,x3,x4}, D, {u1,u2,u3});
    if (isnan(u_tet)) {
      u_tet = std::numeric_limits<double>::infinity();
    }
    // face calculations (contains edge update as fallback)
    double u_face1 = tsitsiklis_update_triangle({x1, x2, x4}, D, {u1, u2});
    double u_face2 = tsitsiklis_update_triangle({x1, x3, x4}, D, {u1, u3});
    double u_face3 = tsitsiklis_update_triangle({x2, x3, x4}, D, {u2, u3});
 
    double u_tri = std::min({u_face1, u_face2, u_face3});
    return std::min(u_tet, u_tri);
  }
}
 
template <int MESH_SIZE>
double LocalSolver<MESH_SIZE>::tsitsiklis_update_line(const std::array<SF::Point, 2>& base,
                                                      const SF::dmat<double>&         D,
                                                      const double&                   value)
{
  // Compute the difference vector: a1 = x2 - x1
  SF::Point a1 = base[1] - base[0];
  // Return updated value at x2: u1 + norm
  return value + std::sqrt(SF::inner_prod(a1, D * a1));
}
 
template <int MESH_SIZE>
double LocalSolver<MESH_SIZE>::tsitsiklis_update_triangle(const std::array<SF::Point, 3>& base,
                                                          const SF::dmat<double>&         D,
                                                          const std::array<double, 2>&    values)
{
  const SF::Point& x1     = base[0];
  const SF::Point& x2     = base[1];
  const SF::Point& x3     = base[2];
  const double&    u1     = values[0];
  const double&    u2     = values[1];
  double           result = std::numeric_limits<double>::infinity();
 
  SF::Point z1 = x1 - x2;
  SF::Point z2 = x2 - x3;
  double    k  = u1 - u2;
 
  // Squared Mahalanobis norms
  const SF::Point Dz1 = D * z1;
  const SF::Point Dz2 = D * z2;
  double p11 = SF::inner_prod(z1, Dz1);
  double p12 = SF::inner_prod(z1, Dz2);
  double p22 = SF::inner_prod(z2, Dz2);
 
  double denominator = p11 - k * k;
  double sqrt_val    = (p11 * p22 - p12 * p12) / denominator;
 
  if (denominator > 1e-12) { // avoid dividing by zero or near zero -> k very small, likely due to collapsed triangle/coinciding points
    const double sqrt_val = (p11 * p22 - p12 * p12) / denominator;
    if (sqrt_val >= 0.0) { // only real solutions are considered
      const double rhs    = k * std::sqrt(sqrt_val);
      double       alpha1 = -(p12 + rhs) / p11;
      double       alpha2 = -(p12 - rhs) / p11;
 
      alpha1 = std::clamp(alpha1, 0.0, 1.0);
      alpha2 = std::clamp(alpha2, 0.0, 1.0);
 
      for (double alpha : {alpha1, alpha2}) {
        SF::Point x_interp = x1 * alpha + x2 * (1.0 - alpha);
        SF::Point dist     = x3 - x_interp;
        double          norm_D   = std::sqrt(SF::inner_prod(dist, D * dist));
        double          u3       = alpha * u1 + (1.0 - alpha) * u2 + norm_D;
        result                   = std::min(result, u3);
      }
    }
  }
 
  // Fallback: point-based update if square root was invalid or denominator is too small
  double u_edge1 = tsitsiklis_update_line({x1,x3}, D, u1);
  double u_edge2 = tsitsiklis_update_line({x2,x3}, D, u2);
 
  return std::min({result, u_edge1, u_edge2});
}
 
template <int MESH_SIZE>
double LocalSolver<MESH_SIZE>::tsitsiklis_update_tetra(const std::array<SF::Point, 4>& base,
                                                       const SF::dmat<double>&         D,
                                                       const std::array<double, 3>&    values)
{
  const SF::Point& x1 = base[0];
  const SF::Point& x2 = base[1];
  const SF::Point& x3 = base[2];
  const SF::Point& x4 = base[3];
 
  const double& u1 = values[0];
  const double& u2 = values[1];
  const double& u3 = values[2];
 
  // edge vectors
  const SF::Point y1 = x3 - x1;
  const SF::Point y2 = x3 - x2;
  const SF::Point y3 = x4 - x3;
 
  const double k1 = u1 - u3;
  const double k2 = u2 - u3;
 
  // Matrix products (squared norms and dot products under metric D)
  const SF::Point Dy1 = D * y1;
  const SF::Point Dy2 = D * y2;
  const SF::Point Dy3 = D * y3;
  const double r11 = SF::inner_prod(y1, Dy1);
  const double r12 = SF::inner_prod(y1, Dy2);
  const double r13 = SF::inner_prod(y1, Dy3);
  const double r21 = r12;
  const double r22 = SF::inner_prod(y2, Dy2);
  const double r23 = SF::inner_prod(y2, Dy3);
  const double r31 = r13;
  const double r32 = r23;
 
  const double A1 = k2 * r11 - k1 * r12;
  const double A2 = k2 * r21 - k1 * r22;
  const double B  = k2 * r31 - k1 * r32;
  const double k  = k1 - (A1 / A2) * k2;
  const SF::Point z1 = y1 - (A1 / A2) * y2;
  const SF::Point z2 = y3 - (B / A2) * y2;
 
  // compute quadratic equation
  const SF::Point Dz1 = D * z1;
  const SF::Point Dz2 = D * z2;
  const double p11 = SF::inner_prod(z1, Dz1);
  const double p12 = SF::inner_prod(z1, Dz2);
  const double p22 = SF::inner_prod(z2, Dz2);
  const double denominator = p11 - k*k;
  const double sqrt_val = (p11 * p22 - (p12 * p12)) / denominator;
  const double rhs = k * std::sqrt(sqrt_val);
 
  double alpha1 = -(p12 + rhs) / p11;
  double alpha2 = -(B + alpha1 * A1) / A2;
 
  // handle degenerate cases
  const double EPS = 1e-16;
  if ((std::abs(A1) < EPS) && (std::abs(A2) < EPS)) {
    alpha1 = (r12 * r23 - r13 * r22) / (r11 * r22 - (r12 * r12));
    alpha2 = (r12 * r13 - r11 * r23) / (r11 * r22 - (r12 * r12));
  } else if ((std::abs(A1) < EPS) && (std::abs(A2) > EPS)) {
    alpha1 = 0;
    alpha2 = -B / A2;
  } else if ((std::abs(A1) > EPS) && (std::abs(A2) < EPS)) {
    alpha1 = -B / A1;
    alpha2 = 0;
  }
 
  double alpha3 = 1 - alpha1 - alpha2;
  const SF::Point dist = x4 - (alpha1 * x1 + alpha2 * x2 + alpha3 * x3);
  // barycentric coordinate should be inside the tetrahedron
  if (alpha1 < -EPS || alpha2 < -EPS || alpha3 < -EPS ||
    alpha1 > 1.0+EPS || alpha2 > 1.0+EPS || alpha3 > 1.0+EPS) {
    return std::numeric_limits<double>::infinity();
  } else {
    return alpha1 * u1 + alpha2 * u2 + alpha3 * u3 + std::sqrt(SF::inner_prod(dist, D * dist));
  }
}
 
}  // namespace opencarp