stimulate.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 _STIMULATE_H
#define _STIMULATE_H

#include "signals.h"
#include "physics_types.h"
#include "fem_utils.h"

namespace opencarp {

// the following must be kept in jive with electrics.prm
#define EOStim              -2
#define Transmembrane_I      0
#define Extracellular_I      1
#define Extracellular_V      2
#define Extracellular_Ground 3
#define Intracellular_I      4
#define Extracellular_V_OL   5
//#define Mono_Extracellular_V 6
#define Transmembrane_I_Grad 7
#define LAT_Triggered        8
#define Vm_clamp             9
#define Prescribed_Phie      10
#define Ignore_Stim          11

// to be replaced by sets
//#define IsVoltageStim(A)   ((A.stimtype==Extracellular_V) || (A.stimtype==Extracellular_V_OL) || (A.stimtype==Vm_clamp))
#define IsExtraV(A)        ((A.stimtype==Extracellular_V) || (A.stimtype==Extracellular_V_OL))
//#define IsVoltageStim_(A)  ((A.type==Extracellular_V) || (A.type==Extracellular_V_OL) || (A.type==Vm_clamp))
#define IsExtraV_(A)       ((A.type==Extracellular_V) || (A.type==Extracellular_V_OL))

#define SELECTED_STIM(A,B) ((A==B) || ((A==Extracellular_V_OL) && (B==Extracellular_V)))

// flags for ignoring stimuli as function of simulation mode
#define IGNORE_NONE  0
#define NO_EXTRA_GND 1
#define NO_EXTRA_V   2
#define NO_EXTRA_I   4

#define STM_IGNORE_BIDOMAIN    (IGNORE_NONE)                            // bidomain
#define STM_IGNORE_MONODOMAIN  (NO_EXTRA_GND | NO_EXTRA_V | NO_EXTRA_I) // monodomain
#define STM_IGNORE_PSEUDO_BIDM (NO_EXTRA_V | NO_EXTRA_I)                // pseudo-bidomain
#define STM_IGNORE_FLOAT_GND   (NO_EXTRA_GND | NO_EXTRA_V)              // bidomain w/floating ground

// define electrode specification method
#define VOL_BASED_ELEC_DEF   0
#define FILE_BASED_ELEC_DEF  1

// define boundary condition specification method
#define VOL_BASED_BC_DEF   0
#define FILE_BASED_BC_DEF  1

enum waveform_t {squarePulse = 0, truncExpPulse, sinePulse, arbPulse, constPulse, unsetPulse};
const std::string wfLabels [] = {"squarePulse", "truncExpPulse", "sinePulse", "arbPulse"};

//enum stim_t {I_transmembrane = 0, I_extra};
enum stim_t {I_tm=0, I_ex=1, V_ex=2, GND_ex=3, I_in=4, V_ex_ol=5, Illum=6, I_tm_grad=7, I_lat=8, Vm_clmp=9, Phie_pre=10, Ignore_stim=11};
enum stim_domain_t {intra_stim = 1, purk_stim = 2, all_stim = 3}; // matches Stimulus.domain from carp.prm

// stimulus info conveniene functions
void init_stim_info(void);
bool is_voltage(stim_t type);    
bool is_current(stim_t type);    
bool is_dbc(stim_t type);        
bool is_extra(stim_t type);      

class stim_pulse
{
  public:
    double      strength  = 0.0;          
    double      duration  = 0.0;          
    waveform_t  wform     = unsetPulse;   

    sig::time_trace  wave;                

    inline void assign(double _strength, double _duration, double _dt, waveform_t _wform)
    {
      strength = _strength;
      duration = _duration;
      wform    = _wform;
      wave     = sig::time_trace(_duration, _dt);
      wave.set_labels(wfLabels[wform]);
    }

    void setup(int id);
};

class stim_protocol
{
  public:
    double start    = 0.0;      
    int    npls     = 0;        
    double pcl      = 0.0;      
    int    timer_id = -1;       
    int    xtrg_id  = -1;       

    void setup(int idx, std::string name);
};

class stim_physics
{
  public:
    stim_t         type;            
    stim_domain_t  domain;          
    std::string    unit;            
    double         scale;           
    bool           total_current;   

    void setup(int idx);
};

class stim_electrode
{
  public:

  enum def_t {file_based, vol_based_tag, vol_based_shape};

  def_t        definition;
  SF::vector<mesh_int_t> vertices;
  SF::vector<SF_real>  scaling;
  std::string  input_filename;

  void setup(int idx);
};


class stimulus
{
  public:
    int            idx  = -1;          
    std::string    name = "unnamed";   
    stim_pulse     pulse;              
    stim_protocol  ptcl;               
    stim_physics   phys;               
    stim_electrode electrode;          

    void setup(int idx);
    void translate(int id);   

    bool value(double & v) const;

    bool is_active() const;
};

class dbc_manager
{
  public:

  struct dbc_data {
    SF::vector<SF_int>* nod  = nullptr;
    sf_vec*         cntr = nullptr;

    ~dbc_data() {
      if(nod)  delete nod;
      if(cntr) delete cntr;
    }
  };

  sf_mat & mat;
  const SF::vector<stimulus> & stimuli;
  std::map<int, dbc_data*> active_dbc;

  dbc_manager(sf_mat & im, const SF::vector<stimulus> & is) :
              mat(im), stimuli(is)
  {
    recompute_dbcs();
  }

  ~dbc_manager()
  {
    clear_active_dbc();
  }

  bool dbc_update();
  void recompute_dbcs();

  void enforce_dbc_lhs();
  void enforce_dbc_rhs(sf_vec & rhs);


  private:
  inline void clear_active_dbc()
  {
    for(auto & d : active_dbc)
      delete d.second;

    active_dbc.clear();
  }
};

void sample_wave_form(stim_pulse& sp,int idx);

}  // namespace opencarp

#endif