electrics_eikonal.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: Apache-2.0
 
#pragma once
 
#ifndef _EIKONAL_H
#define _EIKONAL_H
 
#include <array>
 
#include "SF_globals.h"
#include "electrics.h"
#include "physics_types.h"
#include "sim_utils.h"
#include "stimulate.h"
 
#include "sf_interface.h"
#include "timers.h"
#include "ionics.h"
#include "basics.h"
 
namespace opencarp
{
struct node_stats {
  enum status { in  = 1,
                out = 0 };
  enum reason { none = 0,
                nbn  = 1,
                conv = 2,
                stim = 3 };
  mesh_int_t  idX;                                                   
  mesh_int_t  cycle     = 0;                                         
  mesh_int_t  idXNB     = std::numeric_limits<Int>::min();           
  SF_real     T_A       = std::numeric_limits<double>::quiet_NaN();  
  SF_real     nbn_T_A   = std::numeric_limits<double>::quiet_NaN();  
  SF_real     T_A_      = std::numeric_limits<double>::quiet_NaN();  
  SF_real     T_R       = std::numeric_limits<double>::quiet_NaN();  
  SF_real     D_I       = std::numeric_limits<double>::quiet_NaN();  
  const char* status    = "out";                                     
  const char* reasonIn  = "-";                                       
  const char* reasonOut = "-";                                       
 
  FILE_SPEC logger = NULL;
 
  ~node_stats()
  {
    f_close(logger);
  }
 
  void init_logger(const char* filename);
  void log_stats(double tm, bool cflg);
  void update_status(enum status s, enum reason r);
};
 
class eikonal_solver
{
 public:
  eikonal_solver() : stimuliRef(nullptr) {}
 
  const double inf = std::numeric_limits<double>::infinity();
 
  // data to dump
  node_stats nodeData;
 
  // Has to be kept in line with param_globals::eik_solve
  enum Idiff_t { FOOT  = 0,
                 GAUSS = 1 };
  Idiff_t Idiff_tech = FOOT;
 
  struct diffusion_current {
    Idiff_t model   = FOOT;
    SF_real A_F     = 1;
    SF_real tau_F   = 1;
    SF_real V_th    = 1;
    SF_real alpha_1 = 1;
    SF_real beta_1  = 1;
    SF_real gamma_1 = 1;
    SF_real alpha_2 = 1;
    SF_real beta_2  = 1;
    SF_real gamma_2 = 1;
    SF_real alpha_3 = 1;
    SF_real beta_3  = 1;
    SF_real gamma_3 = 1;
  };
 
  // Vectors
  SF::vector<diffusion_current> diff_cur;
  SF::vector<SF_real>           T_A, T_R, TA_old, D_I;                          // Time of activation, time of repolarization, time of activation previously calculated, diastolic interval, final conduction velocity
  SF::vector<mesh_int_t>        List, num_changes, nReadded2List, stim_status;  // List of active nodes (L), number of changes per node to enter iterate in the list, number of times reentering the list per node, status per node after stimulus
  SF::vector<mesh_int_t>        StimulusPoints;                                 // Nodes that have initial stimulus
  SF::vector<SF_real>           StimulusTimes;                                  // Times when the initial nodes are stimulated
  mesh_int_t                    Index_currStim;  // index of the row in the table of stimuli and points
  SF_real                       actMIN = 0, actMAX = 0;
 
  // mesh related data
  mesh_int_t              MESH_SIZE;  // number of points per element
  bool                    twoFib;     // true if there is sheet information, false if there is only fiber orientation
  size_t                  num_pts;
  SF::vector<mesh_int_t>  e2n_con;     // Connectivity vector with nodes that belong to each element
  SF::vector<mesh_int_t>  elem_start;  // For the i_th element elem_start[j] stores its starting position in the "connect" vector
  std::vector<mesh_int_t> n2n_connect, n2n_dsp;
  SF::vector<mesh_int_t>  e2n_cnt, n2e_cnt, n2e_con, n2e_dsp;
 
  // Element data storage for anisotropy and CV
  std::vector<SF::dmat<double>> S;
  std::vector<double>          CV_L;
 
  // CV restitution parameters
  SF::vector<SF_real> rho_cvrest, kappa_cvrest, theta_cvrest, denom_cvrest;
  
  // export vectors
  sf_vec*             Idiff = nullptr;
  sf_vec*             AT    = nullptr;
  sf_vec*             RT    = nullptr;
 
  // eikonal solver statistics
  eikonal_solver_stats stats;
 
  ~eikonal_solver()
  {
    if (Idiff) delete Idiff;
    if (AT) delete AT;
    if (RT) delete RT;
    if (stimuliRef) delete stimuliRef;
  }
 
  void init();
 
  bool determine_model_to_run(double& time);
 
  void set_stimuli(SF::vector<stimulus>& stimuli) { stimuliRef = &stimuli; }
 
  void clean_list();
 
  void save_eikonal_state(const char* tsav_ext);
 
  void update_repolarization_times(const Ionics& ion);
 
  void cycFIM();
 
  void FIM();
 
  void update_repolarization_times_from_rd(sf_vec& Vmv, sf_vec& Vmv_old, double time);
 
  void compute_diffusion_current(const double& time, sf_vec& vm);
 
  void init_imp_region_properties();
 
 private:
  // stimulus pointer
  SF::vector<stimulus>* stimuliRef = nullptr;
 
  void compute_bc();
 
  void precompute_squared_anisotropy_metric();
 
  // --- Small helpers for active list management ---
  inline void add_to_active(std::vector<mesh_int_t>& active,
                            std::vector<char>&       in_active,
                            mesh_int_t               node)
  {
    if (!in_active[node]) {
      in_active[node] = 1;
      active.push_back(node);
    }
  }
 
  inline void remove_from_active(std::vector<char>& in_active, mesh_int_t node)
  {
    in_active[node] = 0;
  }
 
  void update_Ta_in_active_list();
 
  SF_real compute_H(mesh_int_t& indX, SF_real& tmpTA);
 
  SF_real compute_coherence(mesh_int_t& indX);
 
  SF_real update(mesh_int_t& indX, SF_real CVrest_factor = 1.0, bool isDREAM = false);
 
  template <int N>
  SF_real update_impl(mesh_int_t& indX, SF_real CVrest_factor = 1.0, bool CheckValidity = false);
 
  template <int N>
  bool is_not_valid_update(const std::array<int, N>& nodeIDs, double time);
 
  bool add_node_neighbor_to_list(SF_real& RT, SF_real& oldTA, SF_real& newTA);
 
  void translate_stim_to_eikonal();
 
  void create_node_to_node_graph();
 
  void load_state_file();
};
 
class Eikonal : public Basic_physic
{
 public:
  enum grid_t { intra_grid = 0,
                extra_grid };
  enum linsys_t { elliptic_sys,
                  parabolic_sys };
 
  MaterialType mtype[2];
 
  SF::vector<stimulus> stimuli;
 
  sf_vec* phie_dummy = nullptr;
 
  Ionics ion;
 
  parabolic_solver parab_solver;
 
  eikonal_solver eik_solver;
  // Has to be kept in line with param_globals::eik_solve
  enum eikonal_t { EIKONAL = 0,
                   REp     = 1,
                   REm     = 2,
                   DREAM   = 3 };
  eikonal_t eik_tech = EIKONAL;
 
  LAT_detector lat;
 
  gvec_data gvec;
 
  igb_output_manager output_manager_time;
  igb_output_manager output_manager_cycle;
 
  bool do_output_eikonal = false;
 
  generic_timing_stats IO_stats;
 
  Eikonal() : ion(intra_elec_msh)
  {
    name = "Eikonal";
  }
 
  void initialize();
 
  void destroy();
  void compute_step();
  void output_step();
 
  inline void output_timings()
  {
    // since ionics are included in electrics, we substract ionic timings from electric timings
    compute_time -= ion.compute_time;
    initialize_time -= ion.initialize_time;
 
    // now we call the base class timings output for Electrics
    Basic_physic::output_timings();
    // and then for ionics
    ion.output_timings();
  }
 
  ~Eikonal();
 
  double timer_val(const int timer_id);
 
  std::string timer_unit(const int timer_id);
 
 private:
  void solve_EIKONAL();
  void solve_RE();
  void solve_DREAM();
  void solve_RD();
 
  void setup_stimuli();
  void stimulate_intracellular();
  void clamp_Vm();
 
  void setup_solvers();
 
  void setup_mappings();
 
  void setup_output();
 
  void dump_matrices();
 
  void checkpointing();
 
  void set_elec_tissue_properties(MaterialType* mtype, Eikonal::grid_t g, FILE_SPEC logger);
};
 
template <int MESH_SIZE>
class LocalSolver
{
 private:
  SF::dmat<double>&                 D;
  std::array<SF::Point, MESH_SIZE>& points;
  std::array<double, MESH_SIZE>&    values;
 
  double tsitsiklis_update_line(const std::array<SF::Point, 2>& base,
                                const SF::dmat<double>&         D,
                                const double&                   value);
  double tsitsiklis_update_triangle(const std::array<SF::Point, 3>& base,
                                    const SF::dmat<double>&         D,
                                    const std::array<double, 2>&    values);
  double tsitsiklis_update_tetra(const std::array<SF::Point, 4>& base,
                                 const SF::dmat<double>&         D,
                                 const std::array<double, 3>&    values);
 
 public:
  LocalSolver(SF::dmat<double>& D, std::array<SF::Point, MESH_SIZE>& points, std::array<double, MESH_SIZE>& values)
      : D(D), points(points), values(values) {}
 
  // This function performs the local solver, returning the solution value for the last point with reference to
  // point order of constructor input
  double solve();
};
 
template class LocalSolver<2>;
 
template class LocalSolver<3>;
 
template class LocalSolver<4>;
 
}  // namespace opencarp
 
#endif