electrics_eikonal.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------
 
#ifndef _EIKONAL_H
#define _EIKONAL_H
 
#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) {}
 
  size_t num_pts;
 
  mesh_int_t numPtsxElem;        // number of points per element
  SF_real    time2stop_eikonal;  // time before the fim stops iterating and assess
  SF_real    CV_l;               // Variable to store the conduction velocity in each update
  mesh_int_t Index_currStim;     // index of the row in the table of stimuli and points
 
  const double inf = std::numeric_limits<double>::infinity();  // infinity
 
  // Variables for solver for triangles
  SF_real A, B, C;
  SF_real lam;
  SF_real ae;
  SF_real be;
  SF_real ce;
  SF_real det;
  SF_real Num, Frac;
  SF_real TA, TA_ori;
  SF_real p1, p2;
 
  // data to dump
  node_stats nodeData;
 
  // Variables for solver for tetrahedra
  SF_real c1, c2, c3, c4, c5, c6;
  SF_real a1, s1, s2, s3, a2, s4, b1, b2;
  SF_real R1, R2, R3, R4, R5;
  SF_real asol, bsol;
  SF_real a, b;
  SF_real derivat_a, derivat_b;
 
  double DI;  // diastolic interval
 
  // 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
 
  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
  bool                    twoFib;        // true if there is sheet information, false if there is only fiber orientation
 
  SF::vector<SF_real> fiber_node, sheet_node;                                                          // Fiber and sheet orientation per node
  SF::vector<SF_real> x_coord, y_coord, z_coord;                                                       // x,y,z coordinates of each node
  SF::vector<SF_real> CVnodes, CVnodes_mod, p_cvrest, k_cvrest, j_cvrest, l_cvrest, cv_ar_t, cv_ar_n;  // CV restitution parameters
  sf_vec*             Idiff = nullptr;
  sf_vec*             AT    = nullptr;
  sf_vec*             RT    = nullptr;
 
  std::vector<mesh_int_t> n2n_connect, n2n_dsp;
  SF::vector<mesh_int_t>  e2n_cnt, n2e_cnt, n2e_con, n2e_dsp;
 
  SF::dmat<double> D, iD;
 
  SF_real actMIN = 0, actMAX = 0;
 
  // eikonal solver stats
  eikonal_solver_stats stats;
 
  bool condign1, condign2, condign3, condign4;
 
  ~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_initial_cv();
 
  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 compute_bc();
 
 private:
  // stimulus pointer
  SF::vector<stimulus>* stimuliRef = nullptr;
 
  void bubble_sort(SF::vector<SF_real>& Array, SF::vector<mesh_int_t>& IndexArray, int& length);
 
  void update_Ta_in_active_list();
 
  SF_real compute_H(mesh_int_t& indX, SF_real& T_A_indX);
 
  SF_real compute_coherence(mesh_int_t& indX);
 
  SF_real compute_F(mesh_int_t& indX, SF_real& T_A_indX);
 
  SF_real solve_tet(mesh_int_t& indx3, std::vector<mesh_int_t>& Elem_cur);
 
  SF_real solve_tri(mesh_int_t& indx, std::vector<mesh_int_t>& Elem_cur);
 
  SF_real solve_edges(mesh_int_t& indX, mesh_int_t& indY);
 
  SF_real mult_VtMV(SF::vector<SF_real>& X, SF::dmat<double>& iD, SF::vector<SF_real>& Y);
 
  SF_real compute_AT_at_node(SF_real& pp, mesh_int_t& indxf1, mesh_int_t& indyf1, mesh_int_t& indzf1);
 
  SF_real 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);
 
  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();
 
  void init_diffusion_current();
};
 
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 prepace();
 
  void set_elec_tissue_properties(MaterialType* mtype, Eikonal::grid_t g, FILE_SPEC logger);
};
 
}  // namespace opencarp
 
#endif