stimulate.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: Copyright (c) NumeriCor GmbH
// SPDX-License-Identifier: LicenseRef-APL-1.1
 
#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 stim.crct.type in openCARP.prm
#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 Illumination         6
// #define Transmembrane_I_Grad 7
// #define LAT_Triggered        8
#define Vm_clamp           9
#define Intracellular_V    10
#define Intracellular_V_OL 11
#define Ignore_Stim        12
#define LeadField_I        13
 
#define IsExtraV(A)        ((A.stimtype==Extracellular_V) || (A.stimtype==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=1, sinePulse=2, arbPulse, constPulse, unsetPulse};
const std::string wfLabels [] = {"squarePulse", "truncExpPulse", "sinePulse", "arbPulse", "constPulse", "unset"};
 
//enum stim_t {I_transmembrane = 0, I_extra};
enum stim_t {I_tm=0, I_ex=1, Phi_ex=2, GND_ex=3, I_in=4, Phi_ex_ol=5, Illum=6, I_tm_grad=7, I_lat=8, Vm_clmp=9, Phi_in=10, Phi_in_ol=11, Ignore_stim=12, LF_I=13};
enum stim_domain_t {intra_stim = 1, purk_stim = 2, all_stim = 3}; // matches Stimulus.domain from openCARP.prm
 
// stimulus info conveniene functions
void init_stim_info(void);
bool is_potential(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;           
  std::string    unit;           
  double         scale;          
  bool           total_current;  
 
  void setup(int idx, mesh_t intra_mesh, mesh_t extra_mesh);
};
 
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;
  bool dump_vtx;
 
  void setup(int idx, mesh_t intra_mesh, mesh_t extra_mesh);
};
 
 
class stimulus
{
  public:
    int            idx  = -1;          
    std::string    name = "unnamed";   
 
    stim_pulse     pulse;              
    stim_protocol  ptcl;               
    stim_physics   phys;               
    stim_electrode electrode;          
 
    // in some cases one wants to construct a stimulus that is not associated to the
    // electric meshes. One example is the EMI solver.
    mesh_t associated_intra_mesh = intra_elec_msh;
    mesh_t associated_extra_mesh = extra_elec_msh;
 
    void setup(int idx);
    void translate(int id);   
 
    bool value(double & v) const;
 
    void dump_vtx_file(int idx);
 
    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