openCARP
Doxygen code documentation for the open cardiac electrophysiology simulator openCARP
electrics_eikonal.cc
Go to the documentation of this file.
1 // ----------------------------------------------------------------------------
2 // openCARP is an open cardiac electrophysiology simulator.
3 //
4 // Copyright (C) 2020 openCARP project
5 //
6 // This program is licensed under the openCARP Academic Public License (APL)
7 // v1.0: You can use and redistribute it and/or modify it in non-commercial
8 // academic environments under the terms of APL as published by the openCARP
9 // project v1.0, or (at your option) any later version. Commercial use requires
10 // a commercial license (info@opencarp.org).
11 //
12 // This program is distributed without any warranty; see the openCARP APL for
13 // more details.
14 //
15 // You should have received a copy of the openCARP APL along with this program
16 // and can find it online: http://www.opencarp.org/license
17 // ----------------------------------------------------------------------------
18 
27 #include "electrics_eikonal.h"
28 #include "basics.h"
29 #include "petsc_utils.h"
30 #include "timers.h"
31 #include "stimulate.h"
32 
33 #include "SF_init.h"
34 
35 namespace opencarp
36 {
37 
39 {
40  if (get_size() > 1) {
41  log_msg(NULL, 5, 0, "DREAM/Eikonal physics currently do not support MPI parallelization. Use openMP instead. Aborting!");
42  EXIT(EXIT_FAILURE);
43  }
44 
45  double t1, t2;
46  get_time(t1);
47 
48  set_dir(OUTPUT);
49 
50  // open logger
51  logger = f_open("eikonal.log", param_globals::experiment != 4 ? "w" : "r");
52 
53  // setup mappings between extra and intra grids, algebraic and nodal,
54  // and between PETSc and canonical orderings
55  setup_mappings();
56 
57  eik_tech = static_cast<Eikonal::eikonal_t>(param_globals::dream.solve);
58 
59  // the ionic physics is currently triggered from inside the Electrics to have tighter
60  // control over it.
61  switch (eik_tech) {
62  case EIKONAL: break;
63  default:
64  ion.logger = logger;
65  ion.initialize();
66  }
67 
68  // set up Intracellular tissue
69  set_elec_tissue_properties(mtype, intra_grid, logger);
70  region_mask(intra_elec_msh, mtype[intra_grid].regions, mtype[intra_grid].regionIDs, true, "gregion_i");
71 
72  // add electrics timer for time stepping, add to time stepper tool (TS)
73  double global_time = user_globals::tm_manager->time;
74  timer_idx = user_globals::tm_manager->add_eq_timer(global_time, param_globals::tend, 0,
75  param_globals::dt, 0, "elec::ref_dt", "TS");
76 
77  // electrics stimuli setup
78  setup_stimuli();
79 
80  // set up the linear equation systems. this needs to happen after the stimuli have been
81  // set up, since we need boundary condition info
82  setup_solvers();
83 
84  // the next setup steps require the solvers to be set up, since they use the matrices
85  // generated by those
86 
87  // initialize the LATs detector
88  switch (eik_tech) {
89  case EIKONAL: break; // not available for pure eikonal solve
90  default:
92  }
93 
94  // prepare the electrics output. we skip it if we do post-processing
95  if (param_globals::experiment != EXP_POSTPROCESS)
96  setup_output();
97 
98  if (param_globals::prepacing_bcl > 0)
99  prepace();
100 
101  this->initialize_time += timing(t2, t1);
102 }
103 
104 void Eikonal::set_elec_tissue_properties(MaterialType* mtype, Eikonal::grid_t g, FILE_SPEC logger)
105 {
106  MaterialType* m = mtype + g;
107 
108  // initialize random conductivity fluctuation structure with PrM values
109  m->regions.resize(param_globals::num_gregions);
110 
111  const char* grid_name = g == Eikonal::intra_grid ? "intracellular" : "extracellular";
112  log_msg(logger, 0, 0, "Setting up %s tissue poperties for %d regions ..", grid_name,
113  param_globals::num_gregions);
114 
115  char buf[64];
116  RegionSpecs* reg = m->regions.data();
117 
118  for (size_t i = 0; i < m->regions.size(); i++, reg++) {
119  if (!strcmp(param_globals::gregion[i].name, "")) {
120  snprintf(buf, sizeof buf, ", gregion_%d", int(i));
121  param_globals::gregion[i].name = dupstr(buf);
122  }
123 
124  reg->regname = strdup(param_globals::gregion[i].name);
125  reg->regID = i;
126  reg->nsubregs = param_globals::gregion[i].num_IDs;
127  if (!reg->nsubregs)
128  reg->subregtags = NULL;
129  else {
130  reg->subregtags = new int[reg->nsubregs];
131 
132  for (int j = 0; j < reg->nsubregs; j++)
133  reg->subregtags[j] = param_globals::gregion[i].ID[j];
134  }
135 
136  // describe material in given region
137  elecMaterial* emat = new elecMaterial();
138  emat->material_type = ElecMat;
139 
140  emat->InVal[0] = param_globals::gregion[i].g_il;
141  emat->InVal[1] = param_globals::gregion[i].g_it;
142  emat->InVal[2] = param_globals::gregion[i].g_in;
143 
144  emat->ExVal[0] = param_globals::gregion[i].g_el;
145  emat->ExVal[1] = param_globals::gregion[i].g_et;
146  emat->ExVal[2] = param_globals::gregion[i].g_en;
147 
148  emat->BathVal[0] = param_globals::gregion[i].g_bath;
149  emat->BathVal[1] = param_globals::gregion[i].g_bath;
150  emat->BathVal[2] = param_globals::gregion[i].g_bath;
151 
152  // convert units from S/m -> mS/um
153  for (int j = 0; j < 3; j++) {
154  emat->InVal[j] *= 1e-3 * param_globals::gregion[i].g_mult;
155  emat->ExVal[j] *= 1e-3 * param_globals::gregion[i].g_mult;
156  emat->BathVal[j] *= 1e-3 * param_globals::gregion[i].g_mult;
157  }
158  reg->material = emat;
159  }
160 
161  if ((g == Eikonal::intra_grid && strlen(param_globals::gi_scale_vec)) ||
162  (g == Eikonal::extra_grid && strlen(param_globals::ge_scale_vec))) {
164  sf_mesh& mesh = get_mesh(mt);
165  const char* file = g == Eikonal::intra_grid ? param_globals::gi_scale_vec : param_globals::ge_scale_vec;
166 
167  size_t num_file_entries = SF::root_count_ascii_lines(file, mesh.comm);
168 
169  if (num_file_entries != mesh.g_numelem)
170  log_msg(0, 4, 0, "%s warning: number of %s conductivity scaling entries does not match number of elements!",
171  __func__, get_mesh_type_name(mt));
172 
173  // set up parallel element vector and read data
174  sf_vec* escale;
175  SF::init_vector(&escale, get_mesh(mt), 1, sf_vec::elemwise);
176  escale->read_ascii(g == Eikonal::intra_grid ? param_globals::gi_scale_vec : param_globals::ge_scale_vec);
177 
178  if (get_size() > 1) {
179  // set up element vector permutation and permute
180  if (get_permutation(mt, ELEM_PETSC_TO_CANONICAL, 1) == NULL) {
182  }
184  sc(*escale, false);
185  }
186 
187  // copy data into SF::vector
188  SF_real* p = escale->ptr();
189  m->el_scale.assign(p, p + escale->lsize());
190  escale->release_ptr(p);
191  }
192 }
193 
194 void Eikonal::setup_mappings()
195 {
196  bool intra_exits = mesh_is_registered(intra_elec_msh), extra_exists = mesh_is_registered(extra_elec_msh);
197  assert(intra_exits);
198  const int dpn = 1;
199 
200  // It may be that another physic (e.g. ionic models) has already computed the intracellular mappings,
201  // thus we first test their existence
202  if (get_scattering(intra_elec_msh, ALG_TO_NODAL, dpn) == NULL) {
203  log_msg(logger, 0, 0, "%s: Setting up intracellular algebraic-to-nodal scattering.", __func__);
205  }
207  log_msg(logger, 0, 0, "%s: Setting up intracellular PETSc to canonical permutation.", __func__);
209  }
210 }
211 
213 {
214  double t1, t2;
215  get_time(t1);
216 
217  switch (eik_tech) {
218  case EIKONAL: solve_EIKONAL(); break;
219  case DREAM: solve_DREAM(); break;
220  default: solve_RE(); break;
221  }
222 
224  // output lin solver stats
226  }
227  this->compute_time += timing(t2, t1);
228 
229  // since the traces have their own timing, we check for trace dumps in the compute step loop
232 }
233 
235 {
236  double t1, t2;
237  get_time(t1);
238 
239  switch (eik_tech) {
240  case EIKONAL: break;
241  default:
243  output_manager_time.write_data(); // does not exist in EIKONAL
244  }
245 
246  if (do_output_eikonal) {
248  do_output_eikonal = false;
249  }
250 
251  double curtime = timing(t2, t1);
252  this->output_time += curtime;
253 
254  IO_stats.calls++;
255  IO_stats.tot_time += curtime;
256 
259 }
260 
265 {
266  // close logger
267  f_close(logger);
268 
269  switch (eik_tech) {
270  case EIKONAL:
272  break;
273  default:
274  // output LAT data
278  ion.destroy();
279  }
280 }
281 
282 void Eikonal::setup_stimuli()
283 {
284  // initialize basic stim info data (used units, supported types, etc)
285  init_stim_info();
286  bool dumpTrace = true;
287 
288  if (dumpTrace) set_dir(OUTPUT);
289 
290  stimuli.resize(param_globals::num_stim);
291  for (int i = 0; i < param_globals::num_stim; i++) {
292  // construct new stimulus
293  stimulus& s = stimuli[i];
294 
295  if (param_globals::stim[i].crct.type != 0 && param_globals::stim[i].crct.type != 9) {
296  // In the Eikonal class, only the intracellular domain is registered.
297  // Therefore, only I_tm and Vm_clmp are compatible.
298  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);
299  EXIT(EXIT_FAILURE);
300  }
301 
303  s.translate(i);
304 
305  s.setup(i);
306 
307  log_msg(NULL, 2, 0, "Only geometry, start time, npls, and bcl of stim[%i] are used", i);
308 
309  if (dumpTrace && get_rank() == 0)
310  s.pulse.wave.write_trace(s.name + ".trc");
311  }
312 
314 }
315 
316 void Eikonal::stimulate_intracellular()
317 {
319 
320  // iterate over stimuli
321  for (stimulus& s : stimuli) {
322  if (s.is_active()) {
323  // for active stimuli, deal with the stimuli-type specific stimulus application
324  switch (s.phys.type) {
325  case I_tm: {
326  if (param_globals::operator_splitting) {
327  apply_stim_to_vector(s, *ps.Vmv, true);
328  } else {
329  SF_real Cm = 1.0;
331  SF_real sc = tm.time_step / Cm;
332 
333  ps.Irhs->set(0.0);
334  apply_stim_to_vector(s, *ps.Irhs, true);
335 
336  *ps.tmp_i1 = *ps.IIon;
337  *ps.tmp_i1 -= *ps.Irhs;
338  *ps.tmp_i1 *= sc; // tmp_i1 = sc * (IIon - Irhs)
339 
340  // add ionic, transmembrane and intracellular currents to rhs
341  if (param_globals::parab_solve != parabolic_solver::EXPLICIT)
342  ps.mass_i->mult(*ps.tmp_i1, *ps.Irhs);
343  else
344  *ps.Irhs = *ps.tmp_i1;
345  }
346  break;
347  }
348 
349  case Illum: {
350  sf_vec* illum_vec = ion.miif->gdata[limpet::illum];
351 
352  if (illum_vec == NULL) {
353  log_msg(0, 5, 0, "Cannot apply illumination stim: global vector not present!");
354  EXIT(EXIT_FAILURE);
355  } else {
356  apply_stim_to_vector(s, *illum_vec, false);
357  }
358 
359  break;
360  }
361 
362  default: break;
363  }
364  }
365  }
366 }
367 
368 void Eikonal::clamp_Vm()
369 {
370  for (stimulus& s : stimuli) {
371  if (s.phys.type == Vm_clmp && s.is_active())
373  }
374 }
375 
376 void Eikonal::setup_output()
377 {
378  int rank = get_rank();
379  SF::vector<mesh_int_t>* restr_i = NULL;
380  SF::vector<mesh_int_t>* restr_e = NULL;
381  set_dir(OUTPUT);
382 
383  setup_dataout(param_globals::dataout_i, param_globals::dataout_i_vtx, intra_elec_msh,
384  restr_i, param_globals::num_io_nodes > 0);
385 
386  if (param_globals::dataout_i) {
387  switch (eik_tech) {
388  case EIKONAL:
389  output_manager_cycle.register_output(eik_solver.AT, intra_elec_msh, 1, param_globals::dream.output.atfile, "ms", restr_i);
390  break;
391  case DREAM:
392  output_manager_time.register_output(parab_solver.Vmv, intra_elec_msh, 1, param_globals::vofile, "mV", restr_i);
393  output_manager_cycle.register_output(eik_solver.AT, intra_elec_msh, 1, param_globals::dream.output.atfile, "ms", restr_i);
394  output_manager_cycle.register_output(eik_solver.RT, intra_elec_msh, 1, param_globals::dream.output.rtfile, "ms", restr_i);
395  if (strcmp(param_globals::dream.output.idifffile, "") != 0) {
396  output_manager_time.register_output(eik_solver.Idiff, intra_elec_msh, 1, param_globals::dream.output.idifffile, "muA/cm2", restr_i);
397  }
398  break;
399  default:
400  output_manager_time.register_output(parab_solver.Vmv, intra_elec_msh, 1, param_globals::vofile, "mV", restr_i);
401  output_manager_cycle.register_output(eik_solver.AT, intra_elec_msh, 1, param_globals::dream.output.atfile, "ms", restr_i);
402  if (strcmp(param_globals::dream.output.idifffile, "") != 0) {
403  output_manager_time.register_output(eik_solver.Idiff, intra_elec_msh, 1, param_globals::dream.output.idifffile, "muA/cm2", restr_i);
404  }
405  }
406  }
407 
409 
410  if (param_globals::num_trace) {
411  sf_mesh& imesh = get_mesh(intra_elec_msh);
412  open_trace(ion.miif, param_globals::num_trace, param_globals::trace_node, NULL, &imesh);
413  }
414 
415  // initialize generic logger for IO timings per time_dt
416  IO_stats.init_logger("IO_stats.dat");
417 }
418 
419 void Eikonal::dump_matrices()
420 {
421  std::string bsname = param_globals::dump_basename;
422  std::string fn;
423 
424  set_dir(OUTPUT);
425 
426  // dump monodomain matrices
427  if (param_globals::parab_solve == 1) {
428  // using Crank-Nicolson
429  fn = bsname + "_Ki_CN.bin";
430  parab_solver.lhs_parab->write(fn.c_str());
431  }
432  fn = bsname + "_Ki.bin";
433  parab_solver.rhs_parab->write(fn.c_str());
434 
435  fn = bsname + "_Mi.bin";
436  parab_solver.mass_i->write(fn.c_str());
437 }
438 
441 double Eikonal::timer_val(const int timer_id)
442 {
443  // determine
444  int sidx = stimidx_from_timeridx(stimuli, timer_id);
445  double val = 0.0;
446  if (sidx != -1) {
447  stimuli[sidx].value(val);
448  } else
449  val = std::nan("NaN");
450 
451  return val;
452 }
453 
456 std::string Eikonal::timer_unit(const int timer_id)
457 {
458  int sidx = stimidx_from_timeridx(stimuli, timer_id);
459  std::string s_unit;
460 
461  if (sidx != -1)
462  // found a timer-linked stimulus
463  s_unit = stimuli[sidx].pulse.wave.f_unit;
464 
465  return s_unit;
466 }
467 
468 void Eikonal::setup_solvers()
469 {
470  set_dir(OUTPUT);
471 
472  switch (eik_tech) {
473  case EIKONAL:
474  eik_solver.init();
476  break;
477  default:
478  parab_solver.init();
479  parab_solver.rebuild_matrices(mtype, *ion.miif, logger);
480  eik_solver.init();
482  if (param_globals::dump2MatLab) dump_matrices();
483  }
484 }
485 
486 void Eikonal::checkpointing()
487 {
489 
490  // regular user selected state save
491  if (tm.trigger(iotm_chkpt_list)) {
492  char save_fnm[1024];
493  const char* tsav_ext = get_tsav_ext(tm.time);
494 
495  snprintf(save_fnm, sizeof save_fnm, "%s.%s.roe", param_globals::write_statef, tsav_ext);
496 
497  ion.miif->dump_state(save_fnm, tm.time, intra_elec_msh, false, GIT_COMMIT_COUNT);
498  eik_solver.save_eikonal_state(tsav_ext);
499  }
500 
501  // checkpointing based on interval
502  if (tm.trigger(iotm_chkpt_intv)) {
503  char save_fnm[1024];
504  snprintf(save_fnm, sizeof save_fnm, "checkpoint.%.1f.roe", tm.time);
505  ion.miif->dump_state(save_fnm, tm.time, intra_elec_msh, false, GIT_COMMIT_COUNT);
506  }
507 }
508 
509 void Eikonal::prepace()
510 {
511  // Default stimulation parameters for the single cells
512  // To be eventually be replaced by user input or info from bench
513  float stimdur = 1.;
514  float stimstr = 60.;
515 
516  log_msg(NULL, 0, 0, "Using activation times from file %s to distribute prepacing states\n",
517  param_globals::prepacing_lats);
518  log_msg(NULL, 1, 0, "Assuming stimulus strength %f uA/uF with duration %f ms for prepacing\n",
519  stimstr, stimdur);
520 
522  limpet::MULTI_IF* miif = elec->ion.miif;
523 
524  const sf_mesh& mesh = get_mesh(intra_elec_msh);
525  sf_vec* read_lats;
526  SF::init_vector(&read_lats, mesh, 1, sf_vec::algebraic);
527 
528  // read in the global distributed vector of all activation times
529  set_dir(INPUT);
530  size_t numread = read_lats->read_ascii(param_globals::prepacing_lats);
531  if (numread == 0) {
532  log_msg(NULL, 5, 0, "Failed reading required LATs! Skipping prepacing!");
533  return;
534  }
535  set_dir(OUTPUT);
536 
538  assert(sc != NULL);
539 
540  // permute in-place to petsc permutation
541  bool forward = false;
542  (*sc)(*read_lats, forward);
543 
544  // take care of negative LAT values
545  {
546  SF_real* lp = read_lats->ptr();
547  for (int i = 0; i < read_lats->lsize(); i++)
548  if (lp[i] < 0.0) lp[i] = param_globals::tend + 10.0;
549 
550  read_lats->release_ptr(lp);
551  }
552 
553  // make LATs relative and figure out the first LAT
554  // so we know when to save state of each point
555  SF_real LATmin = read_lats->min();
556 
557  if (LATmin < 0.0) {
558  log_msg(0, 3, 0, "LAT data is not complete. Skipping prepacing.");
559  return;
560  }
561 
562  SF_real offset = floor(LATmin / param_globals::prepacing_bcl) * param_globals::prepacing_bcl;
563  SF_real last_tm = param_globals::prepacing_bcl * param_globals::prepacing_beats;
564 
565  // compute read_lats[i] = last_tm - (read_lats[i] - offset)
566  *read_lats += -offset;
567  *read_lats *= -1.;
568  *read_lats += last_tm;
569 
570  miif->getRealData();
571  SF_real* save_tm = read_lats->ptr();
572  SF_real* vm = miif->gdata[limpet::Vm]->ptr();
573 
574  for (int ii = 0; ii < miif->N_IIF; ii++) {
575  if (!miif->N_Nodes[ii]) continue;
576 
577  // create sorted array of save times.
578  SF::vector<SF::mixed_tuple<double, int>> sorted_save(miif->N_Nodes[ii]); // v1 = time, v2 = index
579  for (int kk = 0; kk < miif->N_Nodes[ii]; kk++) {
580  sorted_save[kk].v1 = save_tm[miif->NodeLists[ii][kk]];
581  sorted_save[kk].v2 = kk;
582  }
583  std::sort(sorted_save.begin(), sorted_save.end());
584 
585  size_t lastidx = sorted_save.size() - 1;
586  int paced = sorted_save[lastidx].v2; // IMP index of latest node
587  int csav = 0;
588 
589  for (double t = 0; t < sorted_save[lastidx].v1; t += param_globals::dt) {
590  if (fmod(t, param_globals::prepacing_bcl) < stimdur &&
591  t < param_globals::prepacing_bcl * param_globals::prepacing_beats - 1)
592  miif->ldata[ii][limpet::Vm][paced] += stimstr * param_globals::dt;
593 
594  compute_IIF(*miif->IIF[ii], miif->ldata[ii], paced);
595 
596  // Vm update always happens now outside of the imp
597  miif->ldata[ii][limpet::Vm][paced] -= miif->ldata[ii][limpet::Iion][paced] * param_globals::dt;
598  vm[miif->NodeLists[ii][paced]] = miif->ldata[ii][limpet::Vm][paced];
599 
600  while (csav < miif->N_Nodes[ii] - 1 && t >= sorted_save[csav].v1)
601  dup_IMP_node_state(*miif->IIF[ii], paced, sorted_save[csav++].v2, miif->ldata[ii]);
602  }
603 
604  // get nodes which may be tied for last
605  while (csav < miif->N_Nodes[ii] - 1)
606  dup_IMP_node_state(*miif->IIF[ii], paced, sorted_save[csav++].v2, miif->ldata[ii]);
607  // ipdate global Vm vector
608  for (int k = 0; k < miif->N_Nodes[ii]; k++) vm[miif->NodeLists[ii][k]] = miif->ldata[ii][limpet::Vm][k];
609  }
610 
611  read_lats->release_ptr(save_tm);
612  miif->gdata[limpet::Vm]->release_ptr(vm);
613  miif->releaseRealData();
614 }
615 
616 void Eikonal::solve_EIKONAL()
617 {
618  if (user_globals::tm_manager->time == 0) {
619  eik_solver.FIM();
622  do_output_eikonal = true;
623  }
624 }
625 
626 void Eikonal::solve_RE()
627 {
628  if (user_globals::tm_manager->time == 0) {
629  eik_solver.FIM();
632  do_output_eikonal = true;
633  }
634  solve_RD();
635 }
636 
637 void Eikonal::solve_DREAM()
638 {
640  solve_RD();
641  } else {
642  if (user_globals::tm_manager->time != 0) {
644  }
645 
646  eik_solver.cycFIM();
648 
651  do_output_eikonal = true;
652  }
653 }
654 
655 void Eikonal::solve_RD()
656 {
657  // if requested, we checkpoint the current state
658  checkpointing();
659 
660  // activation checking
661  const double time = user_globals::tm_manager->time,
662  time_step = user_globals::tm_manager->time_step;
663 
664  lat.check_acts(time);
665  lat.check_quiescence(time, time_step);
666 
667  clamp_Vm();
668 
669  *parab_solver.old_vm = *parab_solver.Vmv; // needed in step D of DREAM
670 
671  // compute ionics update
672  ion.compute_step();
673 
674  // Compute the I diff current
675  *eik_solver.Idiff *= 0;
678 
679  if (eik_tech == REp) {
681  }
682 
683  switch (eik_tech) {
685  default: break;
686  }
687 
688  clamp_Vm();
689 }
690 
692 {
693  if (param_globals::output_level > 1) log_msg(0, 0, 0, "\n *** Initializing Eikonal Solver ***\n");
694 
695  stats.init_logger("eik_stats.dat");
696 
697  if (param_globals::dream.output.debugNode >= 0) {
698  nodeData.idX = param_globals::dream.output.debugNode;
699  char buf[256];
700  snprintf(buf, sizeof buf, "node_%d.dat", nodeData.idX);
701  nodeData.init_logger(buf);
702  }
703 
704  // currently only used for output
705  const sf_mesh& mesh = get_mesh(intra_elec_msh);
706  const int dpn = 1;
707  SF::init_vector(&Idiff, mesh, dpn, sf_vec::algebraic);
708  SF::init_vector(&AT, mesh, dpn, sf_vec::algebraic);
709  SF::init_vector(&RT, mesh, dpn, sf_vec::algebraic);
710 
711  num_pts = mesh.l_numpts;
712 
713  e2n_con = mesh.con; // Connectivity vector with nodes that belong to each element
714  // The 4 nodes from index 4j to 4j+3 belong to the same element for 0<=j<num elements
715  elem_start = mesh.dsp; // For the i_th element elem_start[j] stores its starting position in the "connect" vector
716 
717  cnt_from_dsp(elem_start, e2n_cnt);
718  transpose_connectivity(e2n_cnt, e2n_con, n2e_cnt, n2e_con);
719  dsp_from_cnt(n2e_cnt, n2e_dsp);
720 
721  x_coord.resize(num_pts);
722  y_coord.resize(num_pts);
723  z_coord.resize(num_pts);
724  fiber_node.resize(num_pts * 3, 0);
725  sheet_node.resize(num_pts * 3, 0);
726 
727  for (int i = 0; i < ((mesh.xyz.size() + 1) / 3); i++) {
728  x_coord[i] = mesh.xyz[3 * i + 0];
729  y_coord[i] = mesh.xyz[3 * i + 1];
730  z_coord[i] = mesh.xyz[3 * i + 2];
731  }
732 
733  diff_cur.resize(num_pts);
734 
735  D.set_size(3, 3);
736  iD.set_size(3, 3);
737 
738  List.resize(num_pts, 0); // List of active nodes
739  T_A.resize(num_pts, inf); // Activation times
740  D_I.resize(num_pts, 0); // Diastolic Intervals
741  T_R.resize(num_pts, -400); // Recovery times
742  TA_old.resize(num_pts, inf); // Activation Time is the previous cycle
743  num_changes.resize(num_pts, 0); // Number of Times entered in the List per current LAT
744  nReadded2List.resize(num_pts, 0); // Number of Times entered in the List per current LAT
745  stim_status.resize(num_pts, 0); // Status of nodes
746 
747  CVnodes.resize(num_pts, 0);
748  p_cvrest.resize(num_pts);
749  k_cvrest.resize(num_pts);
750  j_cvrest.resize(num_pts);
751  l_cvrest.resize(num_pts);
752  cv_ar_t.resize(num_pts, -1);
753  cv_ar_n.resize(num_pts, -1);
754 
755  StimulusPoints.resize(num_pts * 10, -1); // Nodes that have initial stimulus
756  StimulusTimes.resize(num_pts * 10, -1); // Times when the initial nodes are stimulated
757 
758  switch (mesh.type[0]) {
759  case 0:
760  numPtsxElem = 4;
761  break;
762 
763  case 6:
764  numPtsxElem = 3;
765  break;
766 
767  case 7:
768  numPtsxElem = 2;
769  break;
770 
771  default:
772  log_msg(0, 5, 0, "Error: Type of element is not compatible with this version of the eikonal model. Use tetrahedra or triangles");
773  EXIT(EXIT_FAILURE);
774  break;
775  }
776 
777  if (strlen(param_globals::start_statef) > 0) load_state_file();
778 
779  create_node_to_node_graph();
780 
781  translate_stim_to_eikonal();
782 
783  init_diffusion_current();
784 }
785 
786 void eikonal_solver::bubble_sort(SF::vector<SF_real>& Array, SF::vector<mesh_int_t>& IndexArray, int& length)
787 {
788  SF_real Temp_a;
789  mesh_int_t Temp_i;
790  for (int i = 0; i < length - 1; i++) {
791  for (int j = 0; j < length - i - 1; j++) {
792  if (Array[j] > Array[j + 1]) {
793  Temp_a = Array[j];
794  Array[j] = Array[j + 1];
795  Array[j + 1] = Temp_a;
796  Temp_i = IndexArray[j];
797  IndexArray[j] = IndexArray[j + 1];
798  IndexArray[j + 1] = Temp_i;
799  }
800  }
801  }
802 }
803 
805 {
806  bool act_is_in_safety_window, empty_list_and_illegal_stimulus, empty_stimulus, stimulus_is_in_safety_window;
807 
808  act_is_in_safety_window = (actMIN > param_globals::dream.tau_s) && (actMIN - time > param_globals::dream.tau_s);
809  empty_stimulus = (StimulusPoints[Index_currStim] == -1);
810  stimulus_is_in_safety_window = (StimulusTimes[Index_currStim] - time > param_globals::dream.tau_s);
811  empty_list_and_illegal_stimulus = (sum(List) == 0) && (empty_stimulus || stimulus_is_in_safety_window);
812 
813  return act_is_in_safety_window || empty_list_and_illegal_stimulus;
814 }
815 
817 {
818  const sf_mesh& mesh = get_mesh(intra_elec_msh);
819  SF::vector<mesh_int_t> vertices;
820  for (int num_g = 0; num_g < param_globals::num_gregions; num_g++) {
821  for (int num_id = 0; num_id <= param_globals::gregion[num_g].num_IDs; num_id++) {
822  indices_from_region_tag(vertices, mesh, param_globals::gregion[num_g].ID[num_id]);
823  for (int i = 0; i < vertices.size(); i++) {
824  CVnodes[vertices[i]] = param_globals::gregion[num_g].dream.vel_l;
825  cv_ar_t[vertices[i]] = param_globals::gregion[num_g].dream.vel_l / param_globals::gregion[num_g].dream.vel_t;
826  if (twoFib) {
827  cv_ar_n[vertices[i]] = param_globals::gregion[num_g].dream.vel_l / param_globals::gregion[num_g].dream.vel_n;
828  }
829  }
830  }
831  }
832  CVnodes_mod = CVnodes;
833 
834  for (int num_i = 0; num_i < param_globals::num_imp_regions; num_i++) {
835  for (int num_id = 0; num_id <= param_globals::imp_region[num_i].num_IDs; num_id++) {
836  indices_from_region_tag(vertices, mesh, param_globals::imp_region[num_i].ID[num_id]);
837  for (int i = 0; i < vertices.size(); i++) {
838  p_cvrest[vertices[i]] = param_globals::imp_region[num_i].dream.CVrest.rho;
839  k_cvrest[vertices[i]] = param_globals::imp_region[num_i].dream.CVrest.kappa;
840  j_cvrest[vertices[i]] = param_globals::imp_region[num_i].dream.CVrest.theta;
841  l_cvrest[vertices[i]] = param_globals::imp_region[num_i].dream.CVrest.psi;
842  }
843  }
844  }
845 
846  if (strlen(param_globals::cv_scale_vec)) {
847  mesh_t mt = intra_elec_msh;
848  const char* file = param_globals::cv_scale_vec;
849 
850  size_t num_file_entries = SF::root_count_ascii_lines(file, mesh.comm);
851 
852  if (num_file_entries != mesh.g_numpts)
853  log_msg(0, 2, 0, "%s warning: number of %s conductivity scaling entries does not match number of points!,%i,%i",
854  __func__, get_mesh_type_name(mt), num_file_entries, mesh.g_numpts);
855 
856  // set up parallel point vector and read data
857  sf_vec* escale;
858  SF::init_vector(&escale, mesh, 1, sf_vec::nodal);
859  escale->read_ascii(file);
860 
861  if (get_size() > 1) {
862  // set up point vector permutation and permute
863  if (get_permutation(mt, PETSC_TO_CANONICAL, 1) == NULL) {
865  }
867  sc(*escale, false);
868  }
869 
870  // copy data into SF::vector
871  SF::vector<double> nd_scale;
872  SF_real* p = escale->ptr();
873  nd_scale.assign(p, p + escale->lsize());
874  escale->release_ptr(p);
875 
876  for (int num_g = 0; num_g < param_globals::num_gregions; num_g++) {
877  for (int num_id = 0; num_id < param_globals::gregion[num_g].num_IDs; num_id++) {
878  indices_from_region_tag(vertices, mesh, param_globals::gregion[num_g].ID[num_id]);
879  for (int i = 0; i < vertices.size(); i++) {
880  CVnodes[vertices[i]] *= nd_scale[vertices[i]];
881  }
882  }
883  }
884  CVnodes_mod = CVnodes;
885  }
886 } // Close set initial CV
887 
888 void eikonal_solver::update_Ta_in_active_list()
889 {
890  bool First_in_List = 1;
891 
892  for (size_t j = 0; j < List.size(); j++) {
893  if (T_A[j] < 0) continue;
894  if (List[j] == 1 && First_in_List) {
895  actMIN = T_A[j];
896  actMAX = T_A[j];
897  First_in_List = 0;
898  continue;
899  }
900  if (List[j] == 1) {
901  if (T_A[j] < actMIN) actMIN = T_A[j];
902  if (T_A[j] > actMAX) actMAX = T_A[j];
903  }
904  }
905 }
906 
908 {
909  for (size_t j = 0; j < List.size(); j++) {
910  if (T_A[j] < user_globals::tm_manager->time) {
911  TA_old[j] = T_A[j];
912  stim_status[j] = 0;
913  nReadded2List[j] = 0;
914 
915  if (List[j] == 1) List[j] = 0;
916  }
917  }
918 }
919 
920 SF_real eikonal_solver::compute_H(mesh_int_t& indX, SF_real& T_A_indX)
921 {
922  float delay = 0;
923  if (stim_status[indX] == 2) delay = 5;
924 
925  if (T_A_indX == inf || (T_R[indX] + delay) == -400 || abs(T_A_indX - (TA_old[indX])) < 2) {
926  return CVnodes[indX];
927  } else if ((T_R[indX]) > T_A_indX) {
928  return CVnodes[indX];
929  } else {
930  float DI_indX = T_A_indX - (T_R[indX] + delay);
931  D_I[indX] = DI_indX;
932 
933  float Factor_DI;
934  float denom = log(p_cvrest[indX]) / l_cvrest[indX];
935  Factor_DI = 1 - p_cvrest[indX] * exp(-(DI_indX + k_cvrest[indX]) * denom);
936 
937  if (Factor_DI > 0) {
938  if (DI_indX < j_cvrest[indX]) {
939  return -CVnodes[indX] * Factor_DI;
940  } else {
941  return CVnodes[indX] * Factor_DI;
942  }
943  } else {
944  return 0;
945  }
946  }
947 }
948 
949 SF_real eikonal_solver::compute_coherence(mesh_int_t& indX)
950 {
951  SF_real p = T_A[indX];
952  int iter_ch = 0;
953  SF_real q = compute_F(indX, p);
954  while (abs(p - q) >= param_globals::dream.fim.tol && iter_ch < param_globals::dream.fim.max_coh) {
955  p = q;
956  q = compute_F(indX, p);
957  iter_ch++;
958  }
959  if (compute_H(indX, q) < 0) return T_A[indX];
960  return q;
961 }
962 
964 {
965  bool EmpList = sum(List) == 0;
966 
967  if (EmpList) {
968  for (size_t j = Index_currStim; j < StimulusTimes.size(); j++) {
969  if ((StimulusPoints[j] == -1) || (StimulusTimes[j] > StimulusTimes[Index_currStim])) {
970  Index_currStim = j;
971  break;
972  }
973 
974  mesh_int_t indNode = StimulusPoints[j];
975  SF_real TimeSt = StimulusTimes[j];
976 
977  bool cond2add = add_node_neighbor_to_list(T_R[indNode], T_A[indNode], TimeSt);
978  if (cond2add && compute_H(indNode, TimeSt) > 0) {
979  if (T_A[indNode] < TimeSt) {
980  CVnodes_mod[indNode] = -CVnodes_mod[indNode];
981  }
982  T_A[indNode] = TimeSt;
983  stim_status[indNode] = 1;
984  stats.bc_status = true;
985 
986  if (param_globals::dream.output.debugNode == indNode) {
987  nodeData.T_A = TimeSt;
988  nodeData.update_status(node_stats::in, node_stats::stim);
989  }
990 
991  for (int ii = n2n_dsp[indNode]; ii < n2n_dsp[indNode + 1]; ii++) {
992  mesh_int_t indXNB = n2n_connect[ii];
993  if (List[indXNB] == 0) {
994  List[indXNB] = 1;
995  if (param_globals::dream.output.debugNode == indXNB) {
996  nodeData.update_status(node_stats::in, node_stats::nbn);
997  nodeData.idXNB = indNode;
998  nodeData.nbn_T_A = TimeSt;
999  }
1000  }
1001  }
1002  }
1003 
1004  if (cond2add && compute_H(indNode, TimeSt) <= 0) {
1005  stim_status[indNode] = 2;
1006  }
1007  }
1008 
1009  // After NB of initial points are assigned remove initial points from the list if they were added in the previous loop.
1010  for (size_t j = 0; j < StimulusPoints.size(); j++) {
1011  if ((StimulusPoints[j] == -1) && (StimulusTimes[j] > StimulusTimes[Index_currStim])) continue;
1012  if (List[StimulusPoints[j]] == 1) {
1013  List[StimulusPoints[j]] = 0;
1014  if (param_globals::dream.output.debugNode == StimulusPoints[j]) nodeData.update_status(node_stats::out, node_stats::stim);
1015  }
1016  }
1017 
1018  if (StimulusTimes[Index_currStim] != -1) actMAX = StimulusTimes[Index_currStim];
1019  } else {
1020  for (size_t i = Index_currStim; i < StimulusTimes.size(); i++) {
1021  mesh_int_t indNode = StimulusPoints[i];
1022  SF_real TimeSt = StimulusTimes[i];
1023 
1024  if (indNode == -1) continue;
1025 
1026  if (TimeSt <= actMAX) {
1027  bool cond2add = add_node_neighbor_to_list(T_R[indNode], T_A[indNode], TimeSt);
1028  if (cond2add && compute_H(indNode, TimeSt) > 0) {
1029  T_A[indNode] = TimeSt;
1030  stim_status[indNode] = 1;
1031  stats.bc_status = true;
1032 
1033  if (param_globals::dream.output.debugNode == indNode) {
1034  nodeData.T_A = TimeSt;
1035  nodeData.update_status(node_stats::in, node_stats::stim);
1036  }
1037 
1038  for (int ii = n2n_dsp[indNode]; ii < n2n_dsp[indNode + 1]; ii++) {
1039  mesh_int_t indXNB = n2n_connect[ii];
1040  if (List[indXNB] == 0) {
1041  List[indXNB] = 1;
1042  if (param_globals::dream.output.debugNode == indXNB) {
1043  nodeData.update_status(node_stats::in, node_stats::nbn);
1044  nodeData.idXNB = indNode;
1045  nodeData.nbn_T_A = TimeSt;
1046  }
1047  }
1048  }
1049 
1050  for (size_t j = 0; j < StimulusPoints.size(); j++) {
1051  if ((StimulusPoints[j] == -1) && (StimulusTimes[j] > StimulusTimes[Index_currStim])) continue;
1052  if (List[StimulusPoints[j]] == 1) {
1053  List[StimulusPoints[j]] = 0;
1054  if (param_globals::dream.output.debugNode == StimulusPoints[j]) nodeData.update_status(node_stats::out, node_stats::stim);
1055  }
1056  }
1057  }
1058  if (cond2add && compute_H(indNode, TimeSt) <= 0) {
1059  stim_status[indNode] = 2;
1060  }
1061 
1062  } else {
1063  Index_currStim = i;
1064 
1065  break;
1066  }
1067  }
1068  }
1069 
1070 } // end compute stimulus
1071 
1073 {
1074  int niter = 0;
1075  double t1, t0;
1076  get_time(t0);
1077 
1078  for (size_t i = Index_currStim; i < StimulusTimes.size(); i++) {
1079  mesh_int_t indNode = StimulusPoints[i];
1080  SF_real TimeSt = StimulusTimes[i];
1081  if (StimulusPoints[i] == -1) break;
1082 
1083  T_A[indNode] = TimeSt;
1084 
1085  for (size_t j = n2n_dsp[indNode]; j < n2n_dsp[indNode + 1]; j++) {
1086  if (T_A[n2n_connect[j]] == inf) {
1087  List[n2n_connect[j]] = 1;
1088  }
1089  }
1090 
1091  List[indNode] = 0;
1092  }
1093 
1094  do {
1095  for (mesh_int_t indX = 0; indX < List.size(); indX++) {
1096  if (List[indX] == 0) continue;
1097 
1098  SF_real p = T_A[indX];
1099  SF_real q = compute_F(indX, p);
1100  T_A[indX] = q;
1101 
1102  if (abs(p - q) < param_globals::dream.fim.tol) {
1103  for (int ii = n2n_dsp[indX]; ii < n2n_dsp[indX + 1]; ii++) {
1104  mesh_int_t indXNB = n2n_connect[ii];
1105  if (List[indXNB] == 1) continue;
1106  p = T_A[indXNB];
1107  q = compute_F(indXNB, p);
1108  if (p > q) {
1109  T_A[indXNB] = q;
1110  List[indXNB] = 1;
1111  }
1112  }
1113  List[indX] = 0;
1114  }
1115  }
1116 
1117  update_Ta_in_active_list();
1118  niter++;
1119 
1120  if (niter > param_globals::dream.fim.max_iter) break;
1121 
1122  } while (sum(List) > 0);
1123 
1124  // copy for igb output
1125  double* atc = AT->ptr();
1126  for (mesh_int_t i = 0; i < List.size(); i++) {
1127  if (T_A[i] == inf) {
1128  // for better visualization in meshalyzer
1129  atc[i] = -1;
1130  } else {
1131  atc[i] = T_A[i];
1132  }
1133  }
1134  AT->release_ptr(atc);
1135 
1136  // treat solver statistics
1137  auto dur = timing(t1, t0);
1138  stats.slvtime_A += dur;
1139  stats.update_iter(niter);
1140  stats.minAT = actMIN;
1141  stats.maxAT = actMAX;
1142  stats.activeList = sum(List);
1143  stats.update_cli(user_globals::tm_manager->time, false);
1144 }
1145 
1147 {
1148  int niter = 0;
1149  double t1, t0;
1150  get_time(t0);
1151 
1152  if (sum(List) == 0) {
1153  compute_bc();
1154  }
1155 
1156  time2stop_eikonal = param_globals::dream.tau_inc;
1157  float maxadvance = user_globals::tm_manager->time + param_globals::dream.tau_s + param_globals::dream.tau_inc + param_globals::dream.tau_max;
1158 
1159  do {
1160  if (StimulusTimes[Index_currStim] <= maxadvance) {
1161  compute_bc();
1162  }
1163 
1164  for (mesh_int_t indX = 0; indX < List.size(); indX++) {
1165  SF_real p = T_A[indX];
1166  SF_real q;
1167 
1168  if (List[indX] == 0) continue;
1169 
1170  q = compute_coherence(indX);
1171 
1172  T_A[indX] = q;
1173  num_changes[indX] = num_changes[indX] + 1;
1174 
1175  if (q > maxadvance) continue;
1176 
1177  if (abs(p - q) < param_globals::dream.fim.tol || (num_changes[indX] > param_globals::dream.fim.max_iter) || (q == inf || p == inf)) {
1178  for (int ii = n2n_dsp[indX]; ii < n2n_dsp[indX + 1]; ii++) {
1179  mesh_int_t indXNB = n2n_connect[ii];
1180 
1181  if (List[indXNB] == 1) {
1182  continue;
1183  }
1184 
1185  SF_real pNB = T_A[indXNB];
1186  SF_real qNB;
1187 
1188  qNB = compute_coherence(indXNB);
1189 
1190  bool node_is_valid = add_node_neighbor_to_list(T_R[indXNB], pNB, qNB) && qNB != inf && qNB > user_globals::tm_manager->time;
1191  // This second condition is there to be able to add a node to the list when a new valid activation time
1192  // is found but a reentry would be blocked by the L2 parameter. Normally this condition does not make or brake the
1193  // simulation but would leave individual nodes inactivated, which is not ideal.
1194  bool ignore_L2_if_valid = node_is_valid && qNB > pNB && nReadded2List[indXNB] >= param_globals::dream.fim.max_addpt;
1195 
1196  if (node_is_valid && (nReadded2List[indXNB] < param_globals::dream.fim.max_addpt) || ignore_L2_if_valid) {
1197  if (qNB > pNB) {
1198  nReadded2List[indXNB] = 0;
1199  }
1200  T_A[indXNB] = qNB;
1201  nReadded2List[indXNB]++;
1202  num_changes[indXNB] = 0;
1203  List[indXNB] = 1;
1204  if (param_globals::dream.output.debugNode == indXNB) {
1205  nodeData.update_status(node_stats::in, node_stats::nbn);
1206  nodeData.idXNB = indX;
1207  nodeData.nbn_T_A = q;
1208  }
1209  }
1210  }
1211 
1212  List[indX] = 0;
1213  if (param_globals::dream.output.debugNode == indX) {
1214  nodeData.update_status(node_stats::out, node_stats::conv);
1215  }
1216  }
1217  }
1218 
1219  SF_real actMIN_old = actMIN;
1220  SF_real progress_time;
1221 
1222  update_Ta_in_active_list();
1223 
1224  if (actMIN > actMIN_old) {
1225  progress_time = actMIN - actMIN_old;
1226  } else {
1227  progress_time = 0;
1228  }
1229 
1230  time2stop_eikonal -= progress_time;
1231  niter++;
1232 
1233  } while (time2stop_eikonal > 0 && sum(List) > 0);
1234 
1235  if (sum(List) == 0) {
1236  if (StimulusTimes[Index_currStim] != -1) actMIN = StimulusTimes[Index_currStim];
1237  }
1238 
1239  // copy for igb output
1240  double* atc = AT->ptr();
1241  double* rpt = RT->ptr();
1242  for (mesh_int_t i = 0; i < List.size(); i++) {
1243  if (T_A[i] == inf) {
1244  // for better visualization in meshalyzer
1245  atc[i] = -1;
1246  } else {
1247  atc[i] = T_A[i];
1248  }
1249  rpt[i] = T_R[i];
1250  }
1251  AT->release_ptr(atc);
1252  RT->release_ptr(rpt);
1253 
1254  // treat solver statistics
1255  auto dur = timing(t1, t0);
1256  stats.slvtime_A += dur;
1257  stats.update_iter(niter);
1258  stats.minAT = actMIN;
1259  stats.maxAT = actMAX;
1260  stats.activeList = sum(List);
1261  stats.update_cli(user_globals::tm_manager->time, false);
1262 
1263  // treat node stats
1264  if (param_globals::dream.output.debugNode >= 0) {
1265  nodeData.T_A = T_A[nodeData.idX];
1266  nodeData.T_R = T_R[nodeData.idX];
1267  nodeData.D_I = D_I[nodeData.idX];
1268  }
1269 
1270 } // close iterate list
1271 
1272 SF_real eikonal_solver::compute_F(mesh_int_t& indX, SF_real& T_A_indX)
1273 {
1274  const double time = user_globals::tm_manager->time;
1275  TA = inf;
1276  SF::vector<mesh_int_t> Winner_Element(3, -1);
1277  CVnodes_mod[indX] = compute_H(indX, T_A_indX);
1278 
1279  CV_l = abs(CVnodes_mod[indX]);
1280 
1281  if (CV_l == 0) {
1282  return TA;
1283  }
1284 
1285  SF::Point fib, she;
1286 
1287  if (twoFib) {
1288  fib.x = fiber_node[3 * indX + 0];
1289  fib.y = fiber_node[3 * indX + 1];
1290  fib.z = fiber_node[3 * indX + 2];
1291 
1292  she.x = sheet_node[3 * indX + 0];
1293  she.y = sheet_node[3 * indX + 1];
1294  she.z = sheet_node[3 * indX + 2];
1295 
1296  SF_real ani_ratio_t2 = cv_ar_t[indX] * cv_ar_t[indX];
1297  SF_real ani_ratio_n2 = cv_ar_n[indX] * cv_ar_n[indX];
1298  SF::Point n = cross(fib, she);
1299  SF::outer_prod(fib, fib, 1, D[0], false);
1300  SF::outer_prod(she, she, 1 / ani_ratio_t2, D[0], true);
1301  SF::outer_prod(n, n, 1 / ani_ratio_n2, D[0], true);
1302 
1303  iD = invert_3x3(D);
1304 
1305  } else {
1306  fib.x = fiber_node[3 * indX + 0];
1307  fib.y = fiber_node[3 * indX + 1];
1308  fib.z = fiber_node[3 * indX + 2];
1309 
1310  SF_real ani_ratio2 = cv_ar_t[indX] * cv_ar_t[indX];
1311  SF::outer_prod(fib, fib, (ani_ratio2 - 1) / ani_ratio2, D[0]);
1312  D[0][0] += (1 / ani_ratio2), D[1][1] += (1 / ani_ratio2), D[2][2] += (1 / ani_ratio2);
1313  iD = invert_3x3(D);
1314  }
1315 
1316  if (numPtsxElem == 4) {
1317  // Solve for Tetrahedron
1318  for (int e_i = n2e_dsp[indX]; e_i < n2e_dsp[indX + 1]; e_i++) {
1319  int Elem_i = n2e_con[e_i];
1320  mesh_int_t indEle = elem_start[Elem_i];
1321 
1322  std::vector<mesh_int_t> Elem_cur = {e2n_con[indEle], e2n_con[indEle + 1], e2n_con[indEle + 2], e2n_con[indEle + 3]};
1323  std::remove(Elem_cur.begin(), Elem_cur.end(), indX); // intentionally used to move indX to the end of Elem_cur but not to be erased
1324 
1325  condign1 = T_A[Elem_cur[0]] == inf || T_A[Elem_cur[1]] == inf || T_A[Elem_cur[2]] == inf || (T_A[Elem_cur[0]] < time) || (T_A[Elem_cur[1]] < time) || (T_A[Elem_cur[2]] < time);
1326  condign2 = (T_A[Elem_cur[0]] < T_R[Elem_cur[0]]) || (T_A[Elem_cur[1]] < T_R[Elem_cur[1]]) || (T_A[Elem_cur[2]] < T_R[Elem_cur[2]]);
1327  condign3 = CV_l == 0;
1328  condign4 = stim_status[indX] == 2 && (stim_status[Elem_cur[0]] == 1 || stim_status[Elem_cur[1]] == 1 || stim_status[Elem_cur[2]] == 1);
1329 
1330  if (!(condign1 || condign2 || condign3 || condign4)) {
1331  SF_real temp3 = solve_tet(indX, Elem_cur);
1332 
1333  if (temp3 < TA && temp3 > T_R[indX] && compute_H(indX, temp3) > 0) {
1334  Winner_Element[0] = Elem_cur[0];
1335  Winner_Element[1] = Elem_cur[1];
1336  Winner_Element[2] = Elem_cur[2];
1337  TA = temp3;
1338  }
1339  }
1340 
1341  std::vector<mesh_int_t> Triang = {-1, -1};
1342  // Solve triangles in Tetrahedra
1343  for (int cases = 0; cases < 3; cases++) {
1344  switch (cases) {
1345  case 0:
1346  Triang[0] = Elem_cur[0];
1347  Triang[1] = Elem_cur[1];
1348  break;
1349 
1350  case 1:
1351  Triang[0] = Elem_cur[0];
1352  Triang[1] = Elem_cur[2];
1353  break;
1354 
1355  case 2:
1356  Triang[0] = Elem_cur[1];
1357  Triang[1] = Elem_cur[2];
1358  break;
1359 
1360  default:
1361  break;
1362  }
1363 
1364  condign1 = T_A[Triang[0]] == inf || T_A[Triang[1]] == inf || (T_A[Triang[0]] < time) || (T_A[Triang[1]] < time);
1365  condign2 = (T_A[Triang[0]] < T_R[Triang[0]]) || (T_A[Triang[1]] < T_R[Triang[1]]); //||abs(T_A[Triang[0]] - TA_old[Triang[0]])<2||abs(T_A[Triang[1]] - TA_old[Triang[1]])<2;
1366  condign3 = CV_l == 0;
1367  condign4 = stim_status[indX] == 2 && (stim_status[Triang[0]] == 1 || stim_status[Triang[1]] == 1);
1368 
1369  if (!(condign1 || condign2 || condign3 || condign4)) {
1370  SF_real temp2 = solve_tri(indX, Triang);
1371 
1372  if (temp2 < TA && temp2 > T_R[indX] && compute_H(indX, temp2) > 0) {
1373  Winner_Element[0] = Triang[0];
1374  Winner_Element[1] = Triang[1];
1375  Winner_Element[2] = -1;
1376  TA = temp2;
1377  }
1378  }
1379  }
1380  }
1381  } // numPtsxElem == 4
1382 
1383  if (numPtsxElem == 3) {
1384  for (int e_i = n2e_dsp[indX]; e_i < n2e_dsp[indX + 1]; e_i++) {
1385  int Elem_i = n2e_con[e_i];
1386  mesh_int_t indEle = elem_start[Elem_i];
1387 
1388  std::vector<mesh_int_t> Elem_cur = {e2n_con[indEle], e2n_con[indEle + 1], e2n_con[indEle + 2]};
1389  std::remove(Elem_cur.begin(), Elem_cur.end(), indX); // intentionally used to move indX to the end of Elem_cur but not to be erased
1390 
1391  condign1 = T_A[Elem_cur[0]] == inf || T_A[Elem_cur[1]] == inf || (T_A[Elem_cur[0]] < time) || (T_A[Elem_cur[1]] < time); //||abs(T_A[Elem_cur[0]] - TA_old[Elem_cur[0]])<2||abs(T_A[Elem_cur[1]] - TA_old[Elem_cur[1]])<2;
1392  condign2 = (T_A[Elem_cur[0]] < T_R[Elem_cur[0]]) || (T_A[Elem_cur[1]] < T_R[Elem_cur[1]]);
1393  condign3 = CV_l == 0;
1394  condign4 = stim_status[indX] == 2 && (stim_status[Elem_cur[0]] == 1 || stim_status[Elem_cur[1]] == 1);
1395 
1396  if (condign1 || condign2 || condign3 || condign4) continue;
1397 
1398  SF_real temp2 = solve_tri(indX, Elem_cur);
1399 
1400  if (TA > temp2 && temp2 > T_R[indX] && compute_H(indX, temp2) > 0) {
1401  TA = temp2;
1402  }
1403  }
1404  } // numPtsxElem == 3
1405 
1406  // Solve for NB Nodes
1407  for (int ii = n2n_dsp[indX]; ii < n2n_dsp[indX + 1]; ii++) {
1408  mesh_int_t indXNB = n2n_connect[ii];
1409  condign4 = stim_status[indX] == 2 && (stim_status[indXNB] == 1);
1410 
1411  if (!(CV_l == 0 || T_A[indXNB] < T_R[indXNB] || T_A[indXNB] == inf || condign4 || T_A[indXNB] < time)) {
1412  SF_real temp1 = solve_edges(indX, indXNB);
1413 
1414  if (temp1 < TA && temp1 > T_R[indX] && compute_H(indX, temp1) > 0) {
1415  Winner_Element[0] = indXNB;
1416  Winner_Element[1] = -1;
1417  Winner_Element[2] = -1;
1418  TA = temp1;
1419  }
1420  }
1421  }
1422 
1423  if (TA == inf) {
1424  Winner_Element[0] = -1;
1425  Winner_Element[1] = -1;
1426  Winner_Element[2] = -1;
1427  TA = T_A[indX];
1428  }
1429 
1430  std::sort(Winner_Element.begin(), Winner_Element.end());
1431  if (compute_H(indX, TA) == 0) TA = inf;
1432 
1433  return TA;
1434 
1435 } // Close Function_F
1436 
1437 SF_real eikonal_solver::solve_tet(mesh_int_t& indx3, std::vector<mesh_int_t>& Elem_cur)
1438 {
1439  SF_real TA_th = inf;
1440 
1441  SF::vector<SF_real> X(3), Y(3), Z(3), W(3);
1442  SF::vector<SF_real> XY(3), XZ(3), XW(3);
1443 
1444  int indy3 = Elem_cur[0], indz3 = Elem_cur[1], indw3 = Elem_cur[2];
1445 
1446  X[0] = x_coord[indx3];
1447  X[1] = y_coord[indx3];
1448  X[2] = z_coord[indx3];
1449 
1450  Y[0] = x_coord[indy3];
1451  Y[1] = y_coord[indy3];
1452  Y[2] = z_coord[indy3];
1453 
1454  Z[0] = x_coord[indz3];
1455  Z[1] = y_coord[indz3];
1456  Z[2] = z_coord[indz3];
1457 
1458  W[0] = x_coord[indw3];
1459  W[1] = y_coord[indw3];
1460  W[2] = z_coord[indw3];
1461 
1462  XY[0] = Y[0] - X[0];
1463  XZ[0] = Z[0] - X[0];
1464  XW[0] = W[0] - X[0];
1465  XY[1] = Y[1] - X[1];
1466  XZ[1] = Z[1] - X[1];
1467  XW[1] = W[1] - X[1];
1468  XY[2] = Y[2] - X[2];
1469  XZ[2] = Z[2] - X[2];
1470  XW[2] = W[2] - X[2];
1471 
1472  c1 = mult_VtMV(XY, iD, XY);
1473  c2 = mult_VtMV(XZ, iD, XZ);
1474  c3 = mult_VtMV(XW, iD, XW);
1475  c4 = 2 * (mult_VtMV(XY, iD, XZ)); //+mult_VtMV(XZ,iD,XY);
1476  c5 = 2 * mult_VtMV(XY, iD, XW); //+mult_VtMV(XW,iD,XY);
1477  c6 = 2 * mult_VtMV(XW, iD, XZ); //+mult_VtMV(XZ,iD,XW);
1478 
1479  a1 = 4 * c2 * c3 - 2 * c2 * c5 - 2 * c3 * c4 + c4 * c6 + c5 * c6 - c6 * c6;
1480  s1 = 4 * c1 * c2 + 4 * c1 * c3 + 4 * c2 * c3 - 4 * c1 * c6 - 4 * c2 * c5 - 4 * c3 * c4 + 2 * c4 * c5 + 2 * c4 * c6 + 2 * c5 * c6 - c4 * c4 - c5 * c5 - c6 * c6;
1481  SF_real s_1 = 1 / s1;
1482  s2 = c1 * c6 * c6 + c2 * c5 * c5 + c3 * c4 * c4 - 4 * c1 * c2 * c3 - c4 * c5 * c6;
1483  s3 = 4 * c1 * c6 - 4 * c1 * c3 - 4 * c2 * c3 - 4 * c1 * c2 + 4 * c2 * c5 + 4 * c3 * c4 - 2 * c4 * c5 - 2 * c4 * c6 - 2 * c5 * c6 + c4 * c4 + c5 * c5 + c6 * c6 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indy3] * c2 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indy3] * c3 - 4 * CV_l * CV_l * T_A[indy3] * T_A[indy3] * c6 + 4 * CV_l * CV_l * T_A[indz3] * T_A[indz3] * c1 + 4 * CV_l * CV_l * T_A[indz3] * T_A[indz3] * c3 - 4 * CV_l * CV_l * T_A[indz3] * T_A[indz3] * c5 + 4 * CV_l * CV_l * T_A[indw3] * T_A[indw3] * c1 + 4 * CV_l * CV_l * T_A[indw3] * T_A[indw3] * c2 - 4 * CV_l * CV_l * T_A[indw3] * T_A[indw3] * c4 - 8 * CV_l * CV_l * T_A[indy3] * T_A[indz3] * c3 - 4 * CV_l * CV_l * T_A[indy3] * T_A[indz3] * c4 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indz3] * c5 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indz3] * c6 - 8 * CV_l * CV_l * T_A[indy3] * T_A[indw3] * c2 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indw3] * c4 - 4 * CV_l * CV_l * T_A[indy3] * T_A[indw3] * c5 + 4 * CV_l * CV_l * T_A[indy3] * T_A[indw3] * c6 - 8 * CV_l * CV_l * T_A[indz3] * T_A[indw3] * c1 + 4 * CV_l * CV_l * T_A[indz3] * T_A[indw3] * c4 + 4 * CV_l * CV_l * T_A[indz3] * T_A[indw3] * c5 - 4 * CV_l * CV_l * T_A[indz3] * T_A[indw3] * c6;
1484  a2 = 2 * CV_l * T_A[indy3] * c2 + 2 * CV_l * T_A[indy3] * c3 - 2 * CV_l * T_A[indy3] * c6 - 2 * CV_l * T_A[indz3] * c3 - CV_l * T_A[indz3] * c4 + CV_l * T_A[indz3] * c5 + CV_l * T_A[indz3] * c6 - 2 * CV_l * T_A[indw3] * c2 + CV_l * T_A[indw3] * c4 - CV_l * T_A[indw3] * c5 + CV_l * T_A[indw3] * c6;
1485  s4 = 4 * c1 * c2 + 4 * c1 * c3 + 4 * c2 * c3 - 4 * c1 * c6 - 4 * c2 * c5 - 4 * c3 * c4 + 2 * c4 * c5 + 2 * c4 * c6 + 2 * c5 * c6 - c4 * c4 - c5 * c5 - c6 * c6;
1486  b1 = 4 * c1 * c3 - 2 * c1 * c6 - 2 * c3 * c4 + c4 * c5 + c5 * c6 - c5 * c5;
1487  b2 = CV_l * T_A[indy3] * c5 - CV_l * T_A[indy3] * c4 - 2 * CV_l * T_A[indy3] * c3 + CV_l * T_A[indy3] * c6 + 2 * CV_l * T_A[indz3] * c1 + 2 * CV_l * T_A[indz3] * c3 - 2 * CV_l * T_A[indz3] * c5 - 2 * CV_l * T_A[indw3] * c1 + CV_l * T_A[indw3] * c4 + CV_l * T_A[indw3] * c5 - CV_l * T_A[indw3] * c6;
1488 
1489  SF_real s2_s3 = s2 / s3;
1490  SF_real LTas = a1 * s_1;
1491  SF_real LTbs = b1 * s_1;
1492 
1493  SF_real RTa = (2 * sqrt(s2_s3) * (a2)) / (s4);
1494  SF_real RTb = (2 * sqrt(s2_s3) * (b2)) / (s4);
1495 
1496  SF_real tol = 1e-24;
1497 
1498  if (abs(s1) > tol && abs(s3) > tol && abs(s4) > tol && (s2 / s3) > 0) {
1499  a = LTas - RTa;
1500  b = LTbs - RTb;
1501 
1502  if ((a > 0) && (a < 1) && (b > 0) && (b < 1) && ((1 - a - b) > 0) && ((1 - a - b) < 1)) {
1503  TA_th = std::min(TA_th, compute_AT_at_node_in_tet(a, b, indx3, indy3, indz3, indw3));
1504  }
1505 
1506  a = LTas + RTa;
1507  b = LTbs + RTb;
1508 
1509  if ((a > 0) && (a < 1) && (b > 0) && (b < 1) && ((1 - a - b) > 0) && ((1 - a - b) < 1)) {
1510  TA_th = std::min(TA_th, compute_AT_at_node_in_tet(a, b, indx3, indy3, indz3, indw3));
1511  }
1512  }
1513 
1514  return TA_th;
1515 }
1516 
1517 SF_real eikonal_solver::solve_tri(mesh_int_t& indx, std::vector<mesh_int_t>& Elem_cur)
1518 {
1519  SF_real TA_tr = inf;
1520 
1521  SF::vector<SF_real> X(3), Y(3), Z(3);
1522  SF::vector<SF_real> XY(3), XZ(3);
1523 
1524  int indy = Elem_cur[0], indz = Elem_cur[1];
1525 
1526  X[0] = x_coord[indx];
1527  X[1] = y_coord[indx];
1528  X[2] = z_coord[indx];
1529 
1530  Y[0] = x_coord[indy];
1531  Y[1] = y_coord[indy];
1532  Y[2] = z_coord[indy];
1533 
1534  Z[0] = x_coord[indz];
1535  Z[1] = y_coord[indz];
1536  Z[2] = z_coord[indz];
1537 
1538  XY[0] = Y[0] - X[0];
1539  XZ[0] = Z[0] - X[0];
1540  XY[1] = Y[1] - X[1];
1541  XZ[1] = Z[1] - X[1];
1542  XY[2] = Y[2] - X[2];
1543  XZ[2] = Z[2] - X[2];
1544 
1545  C = mult_VtMV(XZ, iD, XZ);
1546  A = mult_VtMV(XY, iD, XY) - 2 * mult_VtMV(XZ, iD, XY) + C;
1547  B = 2 * mult_VtMV(XY, iD, XZ) - 2 * C;
1548 
1549  if (abs(T_A[indy] - T_A[indz]) < 1e-5) {
1550  p1 = -B / (2 * A);
1551  if (p1 > 0 && p1 < 1) {
1552  TA_tr = std::min(compute_AT_at_node(p1, indx, indy, indz), TA_tr);
1553  }
1554  }
1555 
1556  lam = ((T_A[indz] - T_A[indy]) * (T_A[indz] - T_A[indy])) * 4 * (CV_l * CV_l);
1557  ae = (4 * (A * A)) - lam * A;
1558  be = 4 * A * B - lam * B;
1559  ce = (B * B) - lam * C;
1560  det = (be * be) - (4 * ae * ce);
1561 
1562  if (det >= 0) {
1563  p1 = (-be + sqrt(det)) / (2 * ae);
1564 
1565  if (p1 > 0 && p1 < 1) {
1566  TA_tr = std::min(compute_AT_at_node(p1, indx, indy, indz), TA_tr);
1567  }
1568 
1569  if (det > 0) {
1570  p2 = (-be - sqrt(det)) / (2 * ae);
1571 
1572  if (p2 > 0 && p2 < 1) {
1573  TA_tr = std::min(compute_AT_at_node(p2, indx, indy, indz), TA_tr);
1574  }
1575  }
1576  }
1577 
1578  return TA_tr;
1579 }
1580 
1581 SF_real eikonal_solver::solve_edges(mesh_int_t& indX, mesh_int_t& indY)
1582 {
1583  SF::vector<SF_real> X(3), Y(3), XY(3);
1584 
1585  X[0] = x_coord[indX];
1586  X[1] = y_coord[indX];
1587  X[2] = z_coord[indX];
1588  Y[0] = x_coord[indY];
1589  Y[1] = y_coord[indY];
1590  Y[2] = z_coord[indY];
1591 
1592  XY[0] = Y[0] - X[0];
1593  XY[1] = Y[1] - X[1];
1594  XY[2] = Y[2] - X[2];
1595 
1596  Num = mult_VtMV(XY, iD, XY);
1597  Num = sqrt(Num);
1598  Frac = Num / CV_l;
1599 
1600  SF_real TAn = T_A[indY] + Frac;
1601 
1602  return TAn;
1603 } // end Solve_Edges
1604 
1605 bool eikonal_solver::add_node_neighbor_to_list(SF_real& RT, SF_real& oldTA, SF_real& newTA)
1606 {
1607  bool isNewTABetweenRTAndOldTA = (RT < newTA) && (newTA < oldTA);
1608  bool isNewTABiggerRTBiggerOldTA = (oldTA < RT) && (RT < newTA);
1609 
1610  return (isNewTABetweenRTAndOldTA) || (isNewTABiggerRTBiggerOldTA);
1611 }
1612 
1613 SF_real eikonal_solver::mult_VtMV(SF::vector<SF_real>& X, SF::dmat<double>& iD, SF::vector<SF_real>& Y)
1614 {
1615  SF_real result = 0;
1616 
1617  if (X[0] == Y[0] && X[1] == Y[1] && X[2] == Y[2]) {
1618  result = iD[0][0] * X[0] * X[0] + 2 * iD[0][1] * X[0] * X[1] + iD[1][1] * X[1] * X[1] + 2 * iD[0][2] * X[0] * X[2] + iD[2][2] * X[2] * X[2] + 2 * iD[1][2] * X[1] * X[2];
1619  } else {
1620  result = iD[0][0] * X[0] * Y[0] + iD[0][1] * (X[0] * Y[1] + X[1] * Y[0]) + iD[1][1] * X[1] * Y[1] + iD[0][2] * (X[0] * Y[2] + X[2] * Y[0]) + iD[2][2] * X[2] * Y[2] + iD[1][2] * (X[1] * X[2] + X[2] * Y[1]);
1621  }
1622 
1623  return result;
1624 
1625 } // mult_VtMV
1626 
1627 SF_real eikonal_solver::compute_AT_at_node(SF_real& pp, mesh_int_t& indxf1, mesh_int_t& indyf1, mesh_int_t& indzf1)
1628 {
1629  if (pp == 0) {
1630  Num = sqrt(C);
1631  Frac = Num / CV_l;
1632  return T_A[indzf1] + Frac;
1633 
1634  } else if (pp == 1) {
1635  Num = sqrt(A + B + C);
1636  Frac = Num / CV_l;
1637  return T_A[indyf1] + Frac;
1638 
1639  } else {
1640  Num = sqrt(A * pp * pp + B * pp + C);
1641  Frac = Num / CV_l;
1642  TA_ori = T_A[indyf1] * pp + T_A[indzf1] * (1 - pp);
1643  return TA_ori + Frac;
1644  }
1645 }
1646 
1647 SF_real eikonal_solver::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)
1648 {
1649  Num = c1 * (pf2 * pf2) + c2 * (qf2 * qf2) + c3 * ((1 - pf2 - qf2) * (1 - pf2 - qf2)) + c4 * (pf2 * qf2) + c5 * (pf2 * (1 - pf2 - qf2)) + c6 * (qf2 * (1 - pf2 - qf2));
1650  Num = sqrt(Num);
1651  Frac = Num / CV_l;
1652  TA_ori = T_A[indyf2] * pf2 + T_A[indzf2] * qf2 + T_A[indwf2] * (1 - pf2 - qf2);
1653  return TA_ori + Frac;
1654 }
1655 
1656 void eikonal_solver::save_eikonal_state(const char* tsav_ext)
1657 {
1658  if (get_rank() == 0) {
1659  FILE* file_writestate;
1660  char buffer_writestate[1024];
1661  snprintf(buffer_writestate, sizeof buffer_writestate, "%s.%s.roe.dat", param_globals::write_statef, tsav_ext);
1662  file_writestate = fopen(buffer_writestate, "w");
1663  for (size_t jjj = 0; jjj < List.size(); jjj++) {
1664  fprintf(file_writestate, "%d %d %d %f %f %f %f \n", List[jjj], num_changes[jjj], nReadded2List[jjj], T_A[jjj], T_R[jjj], TA_old[jjj], D_I[jjj]);
1665  }
1666 
1667  fclose(file_writestate);
1668  }
1669 }
1670 
1672 {
1673  SF_real* old_ptr = Vmv_old.ptr();
1674  SF_real* new_ptr = Vmv.ptr();
1675 
1676  for (int ind_nodes = 0; ind_nodes < num_pts; ind_nodes++) {
1677  double prev_R = T_R[ind_nodes];
1678 
1679  if (old_ptr[ind_nodes] >= param_globals::dream.repol_time_thresh && new_ptr[ind_nodes] < param_globals::dream.repol_time_thresh) {
1680  T_R[ind_nodes] = time;
1681  }
1682  }
1683  Vmv_old.release_ptr(old_ptr);
1684  Vmv.release_ptr(new_ptr);
1685 }
1686 
1688 {
1689  double t1, t0;
1690  get_time(t0);
1691 
1692  #ifdef _OPENMP
1693  int max_threads = omp_get_max_threads();
1694  omp_set_num_threads(1); // with the current setup of this function, it is better to run it in serial.
1695  #endif
1696 
1697  limpet::MULTI_IF* miif = ion.miif;
1698 
1699  const double time = user_globals::tm_manager->time,
1701 
1702  for (int i = 0; i < miif->N_IIF; i++) {
1703  if (!miif->N_Nodes[i]) continue;
1704  limpet::IonIfBase* pIF = ion.miif->IIF[i];
1705  limpet::IonIfBase* IIF_old = pIF->get_type().make_ion_if(pIF->get_target(),
1706  pIF->get_num_node(),
1707  ion.miif->plugtypes[i]);
1708  IIF_old->copy_SVs_from(*pIF, false);
1709 
1710  int current = 0;
1711  do {
1712  int ind_gb = (miif->NodeLists[i][current]);
1713  // Global index of current node;
1714  // save states at the moment
1715  double Vm_old = miif->ldata[i][limpet::Vm][current];
1716  double Iion_old = miif->ldata[i][limpet::Iion][current];
1717  double prev_R = T_R[ind_gb];
1718 
1719  bool cond_2upd = T_R[ind_gb] <= T_A[ind_gb] && T_A[ind_gb] != inf && Vm_old > param_globals::dream.repol_time_thresh;
1720 
1721  if (cond_2upd) {
1722  double elapsed_time = 0;
1723  double prev_Vm = miif->ldata[i][limpet::Vm][current];
1724 
1725  do {
1726  elapsed_time += dt;
1727  pIF->compute(current, current + 1, miif->ldata[i]);
1728  miif->ldata[i][limpet::Vm][current] -= miif->ldata[i][limpet::Iion][current] * param_globals::dt;
1729  if (miif->ldata[i][limpet::Vm][current] < param_globals::dream.repol_time_thresh) {
1730  T_R[ind_gb] = time + elapsed_time;
1731  break;
1732  }
1733  } while (elapsed_time < 600);
1734 
1735  // Restore old states
1736  miif->ldata[i][limpet::Vm][current] = Vm_old;
1737  miif->ldata[i][limpet::Iion][current] = Iion_old;
1738  }
1739  current++;
1740  } while (current < miif->N_Nodes[i]);
1741  pIF->copy_SVs_from(*IIF_old, false);
1742  }
1743 
1744  #ifdef _OPENMP
1745  omp_set_num_threads(max_threads); // restore max threads for parabolic solver
1746  #endif
1747 
1748  // treat solver statistics
1749  auto dur = timing(t1, t0);
1750  stats.slvtime_D += dur;
1751 }
1752 
1754 {
1755  double t1, t0;
1756  get_time(t0);
1757 
1758  SF_real* c = Idiff->ptr();
1759  SF_real* v = vm.ptr();
1760  SF_real dT;
1761  double e_on, e_off;
1762 
1763  int rank = get_rank();
1764  const sf_mesh& mesh = get_mesh(intra_elec_msh);
1765  const SF::vector<mesh_int_t>& alg_nod = mesh.pl.algebraic_nodes();
1766 
1767  for (size_t j = 0; j < alg_nod.size(); j++) {
1768  mesh_int_t loc_nodal_idx = alg_nod[j];
1769  mesh_int_t loc_petsc_idx = local_nodal_to_local_petsc(mesh, rank, loc_nodal_idx);
1770 
1771  if ((time >= T_A[loc_nodal_idx]) && ((T_A[loc_nodal_idx] + 5) >= time) && (List[loc_nodal_idx] == 0)) {
1772  dT = time - T_A[loc_nodal_idx];
1773 
1774  switch (diff_cur[loc_nodal_idx].model) {
1775  case GAUSS:
1776  c[loc_petsc_idx] = diff_cur[loc_nodal_idx].alpha_1 * exp(-((dT - diff_cur[loc_nodal_idx].beta_1) / diff_cur[loc_nodal_idx].gamma_1) * ((dT - diff_cur[loc_nodal_idx].beta_1) / diff_cur[loc_nodal_idx].gamma_1))
1777  + diff_cur[loc_nodal_idx].alpha_2 * exp(-((dT - diff_cur[loc_nodal_idx].beta_2) / diff_cur[loc_nodal_idx].gamma_2) * ((dT - diff_cur[loc_nodal_idx].beta_2) / diff_cur[loc_nodal_idx].gamma_2))
1778  + diff_cur[loc_nodal_idx].alpha_3 * exp(-((dT - diff_cur[loc_nodal_idx].beta_3) / diff_cur[loc_nodal_idx].gamma_3) * ((dT - diff_cur[loc_nodal_idx].beta_3) / diff_cur[loc_nodal_idx].gamma_3));
1779  break;
1780  default:
1781  e_on = (dT >= 0) ? 1.0 : 0.0;
1782  e_off = (v[loc_petsc_idx] < diff_cur[loc_nodal_idx].V_th) ? 1.0 : 0.0;
1783  c[loc_petsc_idx] = diff_cur[loc_nodal_idx].A_F / diff_cur[loc_nodal_idx].tau_F * exp(dT / diff_cur[loc_nodal_idx].tau_F) * e_on * e_off;
1784  }
1785  } // end if
1786  } // end for loop
1787 
1788  Idiff->release_ptr(c);
1789  vm.release_ptr(v);
1790 
1791  // treat solver statistics
1792  auto dur = timing(t1, t0);
1793  stats.slvtime_B += dur;
1794 }
1795 
1796 void eikonal_solver::load_state_file()
1797 {
1798  set_dir(INPUT);
1799  FILE* file_startstate;
1800  char buffer_startstate[strlen(param_globals::start_statef) + 10];
1801 
1802  snprintf(buffer_startstate, sizeof buffer_startstate, "%s.dat", param_globals::start_statef);
1803  file_startstate = fopen(buffer_startstate, "r");
1804 
1805  if (file_startstate == NULL) {
1806  log_msg(NULL, 5, 0, "Not able to open state file: %s", buffer_startstate);
1807  } else if (param_globals::output_level) {
1808  log_msg(NULL, 0, 0, "Open state file for eikonal model: %s", buffer_startstate);
1809  }
1810 
1811  int ListVal, numChangesVal, numChanges2Val;
1812  float T_AVal, T_RVal, TA_oldVal, D_IVal, PCLVal;
1813  int index = 0;
1814 
1815  while (fscanf(file_startstate, "%d %d %d %f %f %f %f", &ListVal, &numChangesVal, &numChanges2Val, &T_AVal, &T_RVal, &TA_oldVal, &D_IVal) == 7) {
1816  List[index] = ListVal;
1817  num_changes[index] = numChangesVal;
1818  nReadded2List[index] = numChanges2Val;
1819  T_A[index] = T_AVal;
1820  T_R[index] = T_RVal;
1821  TA_old[index] = TA_oldVal;
1822  D_I[index] = D_IVal;
1823  ++index;
1824  }
1825  if (param_globals::output_level) log_msg(NULL, 0, 0, "Number of nodes in active list: %i", sum(List));
1826 
1827  fclose(file_startstate);
1828 }
1829 
1830 void eikonal_solver::init_diffusion_current()
1831 {
1832  const sf_mesh& mesh = get_mesh(intra_elec_msh);
1833  SF::vector<mesh_int_t> vertices;
1834 
1835  for (int num_i = 0; num_i < param_globals::num_imp_regions; num_i++) {
1836  for (int num_id = 0; num_id <= param_globals::imp_region[num_i].num_IDs; num_id++) {
1837  indices_from_region_tag(vertices, mesh, param_globals::imp_region[num_i].ID[num_id]);
1838  for (int i = 0; i < vertices.size(); i++) {
1839  diff_cur[vertices[i]].model = static_cast<eikonal_solver::Idiff_t>(param_globals::imp_region[num_i].dream.Idiff.model);
1840  switch (diff_cur[vertices[i]].model) {
1841  case GAUSS:
1842  diff_cur[vertices[i]].alpha_1 = param_globals::imp_region[num_i].dream.Idiff.alpha_i[0];
1843  diff_cur[vertices[i]].alpha_2 = param_globals::imp_region[num_i].dream.Idiff.alpha_i[1];
1844  diff_cur[vertices[i]].alpha_3 = param_globals::imp_region[num_i].dream.Idiff.alpha_i[2];
1845  diff_cur[vertices[i]].beta_1 = param_globals::imp_region[num_i].dream.Idiff.beta_i[0];
1846  diff_cur[vertices[i]].beta_2 = param_globals::imp_region[num_i].dream.Idiff.beta_i[1];
1847  diff_cur[vertices[i]].beta_3 = param_globals::imp_region[num_i].dream.Idiff.beta_i[2];
1848  diff_cur[vertices[i]].gamma_1 = param_globals::imp_region[num_i].dream.Idiff.gamma_i[0];
1849  diff_cur[vertices[i]].gamma_2 = param_globals::imp_region[num_i].dream.Idiff.gamma_i[1];
1850  diff_cur[vertices[i]].gamma_3 = param_globals::imp_region[num_i].dream.Idiff.gamma_i[2];
1851  break;
1852  default:
1853  diff_cur[vertices[i]].A_F = param_globals::imp_region[num_i].dream.Idiff.A_F;
1854  diff_cur[vertices[i]].tau_F = param_globals::imp_region[num_i].dream.Idiff.tau_F;
1855  diff_cur[vertices[i]].V_th = param_globals::imp_region[num_i].dream.Idiff.V_th;
1856  };
1857  }
1858  }
1859  }
1860 }
1861 
1862 void eikonal_solver::translate_stim_to_eikonal()
1863 {
1864  // Set stimulus for eikonal model
1865  if (param_globals::output_level) log_msg(NULL, 0, 0, "Initializing stimuli for eikonal model ...");
1866 
1867  SF::vector<SF_real> StartStimulus;
1868  SF::vector<mesh_int_t> IndexStimulus;
1869  int number_stimulus = param_globals::num_stim;
1870  int total_stimuli = 0;
1871  int count_stm = 0;
1872  int count = 0;
1873  Index_currStim = 0;
1874 
1875  for (stimulus s : *stimuliRef) {
1876  total_stimuli += s.ptcl.npls;
1877  for (int idx_pls = 0; idx_pls < s.ptcl.npls; idx_pls++) {
1878  SF_real start_time = s.ptcl.start + s.ptcl.pcl * idx_pls;
1879  StartStimulus.push_back(start_time);
1880  IndexStimulus.push_back(s.idx);
1881  count_stm++;
1882  }
1883  }
1884  // sort by start times in ascending order
1885  bubble_sort(StartStimulus, IndexStimulus, total_stimuli);
1886 
1887  for (int ind_stim = 0; ind_stim < total_stimuli; ind_stim++) {
1888  stimulus s = (*stimuliRef)[IndexStimulus[ind_stim]];
1889  // Loop over vertices for each defined stim
1890  for (int ind_nodes = 0; ind_nodes < s.electrode.vertices.size(); ind_nodes++) {
1891  StimulusPoints[count] = s.electrode.vertices[ind_nodes];
1892  StimulusTimes[count] = StartStimulus[ind_stim];
1893  count++;
1894  }
1895  }
1896  // Add -1 at the end to indicate we ran out of stimulations
1897  StimulusPoints.resize(count + 2, -1); // Nodes that have initial stimulus
1898  StimulusTimes.resize(count + 2, -1); // Times when the initial nodes are stimulated
1899 }
1900 
1901 void eikonal_solver::create_node_to_node_graph()
1902 {
1903  // Create the Nodes to Nodes Graph (Per Each node the graph finds the NB nodes)
1904  n2n_dsp.resize(num_pts + 1, 0);
1905  const sf_mesh& mesh = get_mesh(intra_elec_msh);
1906  twoFib = mesh.she.size() > 0;
1907 
1908  for (int point_idx = 0; point_idx < num_pts; point_idx++) {
1909  int numNBElem = n2e_dsp[point_idx + 1] - n2e_dsp[point_idx];
1910  std::vector<mesh_int_t> NNNodesxNodes(numNBElem * numPtsxElem);
1911 
1912  for (int eedsp = 0; eedsp < numNBElem; eedsp++) {
1913  int currElem = n2e_con[n2e_dsp[point_idx] + eedsp];
1914 
1915  if (twoFib) {
1916  fiber_node[3 * point_idx] = mesh.fib[3 * currElem];
1917  fiber_node[3 * point_idx + 1] = mesh.fib[3 * currElem + 1];
1918  fiber_node[3 * point_idx + 2] = mesh.fib[3 * currElem + 2];
1919 
1920  sheet_node[3 * point_idx] = mesh.she[3 * currElem];
1921  sheet_node[3 * point_idx + 1] = mesh.she[3 * currElem + 1];
1922  sheet_node[3 * point_idx + 2] = mesh.she[3 * currElem + 2];
1923 
1924  } else {
1925  fiber_node[3 * point_idx] = mesh.fib[3 * currElem];
1926  fiber_node[3 * point_idx + 1] = mesh.fib[3 * currElem + 1];
1927  fiber_node[3 * point_idx + 2] = mesh.fib[3 * currElem + 2];
1928  }
1929 
1930  for (int nndsp = 0; nndsp < numPtsxElem; nndsp++) {
1931  NNNodesxNodes[eedsp * numPtsxElem + nndsp] = e2n_con[elem_start[currElem] + nndsp];
1932  }
1933  }
1934 
1935  std::sort(NNNodesxNodes.begin(), NNNodesxNodes.end());
1936  auto it = std::unique(NNNodesxNodes.begin(), NNNodesxNodes.end());
1937  NNNodesxNodes.erase(it, NNNodesxNodes.end());
1938  it = std::remove(NNNodesxNodes.begin(), NNNodesxNodes.end(), point_idx);
1939  NNNodesxNodes.erase(it, NNNodesxNodes.end());
1940  int max = NNNodesxNodes.size();
1941 
1942  n2n_dsp[point_idx + 1] = n2n_dsp[point_idx] + NNNodesxNodes.size();
1943 
1944  n2n_connect.insert(n2n_connect.end(), NNNodesxNodes.begin(), NNNodesxNodes.end());
1945  }
1946 
1947  n2n_dsp[num_pts] = n2n_connect.size();
1948 }
1949 
1951 {
1952  logger = f_open(filename, "w");
1953 
1954  const char* h1 = " ------ ---------- ---------- ------- ------- | List logic ----- -------- | Neighbor node ---- |";
1955  const char* h2 = " cycle AT old AT RT DI | Status Entry Exit | ID AT |";
1956 
1957  if (logger == NULL)
1958  log_msg(NULL, 3, 0, "%s error: Could not open file %s in %s. Turning off logging.\n",
1959  __func__, filename);
1960  else {
1961  log_msg(logger, 0, 0, "%s", h1);
1962  log_msg(logger, 0, 0, "%s", h2);
1963  }
1964 }
1965 
1966 void node_stats::log_stats(double time, bool cflg)
1967 {
1968  if (!this->logger) return;
1969 
1970  char abuf[256];
1971  char bbuf[256];
1972  char cbuf[256];
1973 
1974  // create nicer output in logger for inf, -inf, nan
1975  std::ostringstream oss_TA, oss_TA_, oss_TR, oss_DI, oss_nbnTA;
1976  oss_TA << this->T_A;
1977  oss_TA_ << this->T_A_;
1978  oss_TR << this->T_R;
1979  oss_DI << this->D_I;
1980  oss_nbnTA << this->nbn_T_A;
1981 
1982  if (this->idXNB == std::numeric_limits<Int>::min()) {
1983  snprintf(cbuf, sizeof cbuf, "%7s %10s", "-", "-");
1984  } else {
1985  snprintf(cbuf, sizeof cbuf, "%7d %10s", this->idXNB, oss_nbnTA.str().c_str());
1986  }
1987 
1988  snprintf(abuf, sizeof abuf, "%6d %10s %10s %7s %7s", this->cycle, oss_TA.str().c_str(), oss_TA_.str().c_str(), oss_TR.str().c_str(), oss_DI.str().c_str());
1989  snprintf(bbuf, sizeof bbuf, "%7s %8s %8s", this->status, this->reasonIn, this->reasonOut);
1990 
1991  unsigned char flag = cflg ? ECHO : 0;
1992  log_msg(this->logger, 0, flag | FLUSH | NONL, "%9.3f %s | %s | %s |\n", time, abuf, bbuf, cbuf);
1993 
1994  this->reasonIn = "-";
1995  this->reasonOut = "-";
1996  this->T_A_ = this->T_A;
1997  this->T_A = std::numeric_limits<double>::quiet_NaN();
1998  this->T_R = std::numeric_limits<double>::quiet_NaN();
1999  this->D_I = std::numeric_limits<double>::quiet_NaN();
2000  this->idXNB = std::numeric_limits<Int>::min();
2001  this->nbn_T_A = std::numeric_limits<double>::quiet_NaN();
2002  this->cycle++;
2003 }
2004 
2006 {
2007  // directly convert enums to strings for logging
2008  const char* stat_str;
2009  const char* reas_str;
2010  switch (s) {
2011  case node_stats::in: stat_str = "in"; break;
2012  case node_stats::out: stat_str = "out"; break;
2013  }
2014 
2015  switch (r) {
2016  case node_stats::none: reas_str = "-"; break;
2017  case node_stats::nbn: reas_str = "nbn"; break;
2018  case node_stats::conv: reas_str = "conv"; break;
2019  case node_stats::stim: reas_str = "stim"; break;
2020  }
2021 
2022  this->status = stat_str;
2023  if (s == in) {
2024  this->reasonIn = reas_str;
2025  } else {
2026  this->reasonOut = reas_str;
2027  }
2028 }
2029 
2030 } // namespace opencarp
constexpr T min(T a, T b)
Definition: ion_type.h:33
SF::scattering * get_permutation(const int mesh_id, const int perm_id, const int dpn)
Get the PETSC to canonical permutation scattering for a given mesh and number of dpn.
int ** NodeLists
local partitioned node lists for each IMP stored
Definition: MULTI_ION_IF.h:201
stim_electrode electrode
electrode geometry
Definition: stimulate.h:168
void init()
Initialize vectors and variables in the eikonal_solver class.
int * subregtags
FEM tags forming this region.
Definition: fem_types.h:97
char * regname
name of region
Definition: fem_types.h:94
void log_stats(double tm, bool cflg)
Definition: timers.cc:93
virtual void release_ptr(S *&p)=0
void FIM()
Standard fast iterative method to solve eikonal equation with active list approach.
double z
Definition: SF_container.h:68
#define FLUSH
Definition: basics.h:311
stim_pulse pulse
stimulus wave form
Definition: stimulate.h:165
char * dupstr(const char *old_str)
Definition: basics.cc:44
description of materal properties in a mesh
Definition: fem_types.h:110
virtual void add_scaled(const abstract_vector< T, S > &vec, S k)=0
physMaterial * material
material parameter description
Definition: fem_types.h:98
sf_vec * Irhs
weighted transmembrane currents
Definition: electrics.h:105
eikonal_solver eik_solver
Solver for the eikonal equation.
LAT_detector lat
the activation time detector
MPI_Comm comm
the parallel mesh is defined on a MPI world
Definition: SF_container.h:392
void set_stimuli(SF::vector< stimulus > &stimuli)
Simple setter for stimulus vector.
igb_output_manager output_manager_cycle
constexpr T max(T a, T b)
Definition: ion_type.h:31
parabolic_solver parab_solver
Solver for the parabolic bidomain equation.
#define PETSC_TO_CANONICAL
Permute algebraic data from PETSC to canonical ordering.
Definition: sf_interface.h:74
const IonType & get_type() const
Gets this IMP&#39;s model type.
Definition: ION_IF.cc:144
physMat_t material_type
ID of physics material.
Definition: fem_types.h:53
vector< T > dsp
connectivity starting index of each element
Definition: SF_container.h:404
sf_vec * phie_dummy
no elliptic solver needed, but we need a dummy for phie to use parabolic solver
overlapping_layout< T > pl
nodal parallel layout
Definition: SF_container.h:417
vector< S > fib
fiber direction
Definition: SF_container.h:408
int calls
# calls for this interval, this is incremented externally
Definition: timers.h:70
virtual void set(const vector< T > &idx, const vector< S > &vals, const bool additive=false, const bool local=false)=0
void update_status(enum status s, enum reason r)
void cnt_from_dsp(const std::vector< T > &dsp, std::vector< T > &cnt)
Compute counts from displacements.
Definition: kdpart.hpp:149
void compute_diffusion_current(const double &time, sf_vec &vm)
Compute the diffusion current approximation I_diff in all activated nodes. Only converged nodes are c...
eikonal_solver_stats stats
void solve(sf_vec &phie_i)
Definition: electrics.cc:1235
void initialize()
Definition: ionics.cc:60
virtual T lsize() const =0
FILE_SPEC logger
The logger of the physic, each physic should have one.
Definition: physics_types.h:64
timer_manager * tm_manager
a manager for the various physics timers
Definition: main.cc:52
double ExVal[3]
extracellular conductivity eigenvalues
Definition: fem_types.h:62
const char * name
The name of the physic, each physic should have one.
Definition: physics_types.h:62
sf_mat * lhs_parab
lhs matrix (CN) to solve parabolic
Definition: electrics.h:111
int add_eq_timer(double istart, double iend, int ntrig, double iintv, double idur, const char *iname, const char *poolname=nullptr)
Add a equidistant step timer to the array of timers.
Definition: timer_utils.cc:78
SF::vector< mesh_int_t > vertices
Definition: stimulate.h:152
size_t read_ascii(FILE *fd)
double time_step
global reference time step
Definition: timer_utils.h:78
void cycFIM()
Implementation of the cyclical fast iterative method used in step A of the DREAM model.
void count(const vector< T > &data, vector< S > &cnt)
Count number of occurrences of indices.
Definition: SF_vector.h:332
const char * get_mesh_type_name(mesh_t t)
get a char* to the name of a mesh type
Definition: sf_interface.cc:46
#define ECHO
Definition: basics.h:308
vector< S > xyz
node cooridnates
Definition: SF_container.h:415
vector< T > con
Definition: SF_container.h:400
Target get_target() const
Definition: ION_IF.h:338
sf_vec * tmp_i1
scratch vector for i-grid
Definition: electrics.h:103
void log_stats(double tm, bool cflg)
mesh_t
The enum identifying the different meshes we might want to load.
Definition: sf_interface.h:58
T & push_back(T val)
Definition: SF_vector.h:283
std::vector< IonTypeList > plugtypes
plugins types for each region
Definition: MULTI_ION_IF.h:215
bool trigger(int ID) const
Definition: timer_utils.h:166
double InVal[3]
intracellular conductivity eigenvalues
Definition: fem_types.h:61
int N_IIF
how many different IIF&#39;s
Definition: MULTI_ION_IF.h:211
void init_stim_info(void)
uses voltage for stimulation
Definition: stimulate.cc:49
double timer_val(const int timer_id)
figure out current value of a signal linked to a given timer
Point and vector struct.
Definition: SF_container.h:65
void close_files_and_cleanup()
close file descriptors
Definition: sim_utils.cc:1527
void init_logger(const char *filename)
Definition: timers.cc:77
void init_vector(SF::abstract_vector< T, S > **vec)
Definition: SF_init.h:99
virtual void copy_SVs_from(IonIfBase &other, bool alloc)=0
Copies the state variables of an IMP.
void region_mask(mesh_t meshspec, SF::vector< RegionSpecs > &regspec, SF::vector< int > &regionIDs, bool mask_elem, const char *reglist)
classify elements/points as belonging to a region
Definition: ionics.cc:362
std::vector< IonIfBase * > IIF
array of IIF&#39;s
Definition: MULTI_ION_IF.h:202
#define ELEM_PETSC_TO_CANONICAL
Permute algebraic element data from PETSC to canonical ordering.
Definition: sf_interface.h:76
bool determine_model_to_run(double &time)
Determine the next model to run in the alternation between RD and Eikonal.
const T * end() const
Pointer to the vector&#39;s end.
Definition: SF_vector.h:128
void initialize()
Initialize the Eikonal class.
generic_timing_stats IO_stats
vec3< V > cross(const vec3< V > &a, const vec3< V > &b)
Definition: vect.h:144
void compute_IIF(limpet::IonIfBase &pIF, limpet::GlobalData_t **impdata, int n)
Definition: ionics.cc:464
void output(vector< int > nodes, IGBheader *h, char *ofname, enum_format of, int t0, int t1, int stride, bool explode, float scale)
Definition: IGBextract.cc:208
SF::scattering * get_scattering(const int from, const int to, const SF::SF_nbr nbr, const int dpn)
Get a scattering from the global scatter registry.
virtual S * ptr()=0
void indices_from_region_tag(SF::vector< mesh_int_t > &idx, const sf_mesh &mesh, const int tag)
Populate vertex data with the vertices of a given tag region.
Definition: fem_utils.cc:155
bool using_legacy_stimuli
flag storing whether legacy stimuli are used
Definition: main.cc:58
void init(sf_vec &vm, sf_vec &phie, int offset, enum physic_t=elec_phys)
initializes all datastructs after electric solver setup
Definition: electrics.cc:1466
virtual void mult(const abstract_vector< T, S > &x, abstract_vector< T, S > &b) const =0
void dump_state(char *, float, opencarp::mesh_t gid, bool, unsigned int)
void transpose_connectivity(const vector< T > &a_cnt, const vector< T > &a_con, vector< T > &b_cnt, vector< T > &b_con)
Transpose CRS matrix graph A into B.
sf_vec * Vmv
global Vm vector
Definition: electrics.h:99
double tot_time
total time, this is incremented externally
Definition: timers.h:72
int mesh_int_t
Definition: SF_container.h:46
void log_stats(double time, bool cflg)
Definition: timers.cc:130
void clean_list()
Clean the list of nodes by resetting their status and tracking changes based on the time step of the ...
virtual void write(const char *filename) const =0
void assemble_sv_gvec(gvec_data &gvecs, limpet::MULTI_IF *miif)
Definition: ionics.cc:637
double time
current time
Definition: timer_utils.h:76
Represents the ionic model and plug-in (IMP) data structure.
Definition: ION_IF.h:138
virtual IonIfBase * make_ion_if(Target target, int num_node, const std::vector< std::reference_wrapper< IonType >> &plugins) const =0
Generate an IonIf object from this type.
void destroy()
Currently we only need to close the file logger.
region based variations of arbitrary material parameters
Definition: fem_types.h:93
GlobalData_t *** ldata
data local to each IMP
Definition: MULTI_ION_IF.h:205
int check_acts(double tm)
check activations at sim time tm
Definition: electrics.cc:1586
double BathVal[3]
bath conductivity eigenvalues
Definition: fem_types.h:63
std::string name
label stimulus
Definition: stimulate.h:164
int timer_idx
the timer index received from the timer manager
Definition: physics_types.h:66
MaterialType mtype[2]
the material types of intra_grid and extra_grid grids.
void write_data()
write registered data to disk
Definition: sim_utils.cc:1480
void compute(int start, int end, GlobalData_t **data)
Perform ionic model computation for 1 time step.
Definition: ION_IF.cc:270
sf_mesh & get_mesh(const mesh_t gt)
Get a mesh by specifying the gridID.
Definition: sf_interface.cc:33
SF::scattering * register_permutation(const int mesh_id, const int perm_id, const int dpn)
Register a permutation between two orderings for a mesh.
void set_initial_cv()
Set the initial conduction velocity (CV) parameters for the nodes in the mesh.
SF::vector< stimulus > stimuli
the electrical stimuli
const char * get_tsav_ext(double time)
Definition: electrics.cc:831
#define ALG_TO_NODAL
Scatter algebraic to nodal.
Definition: sf_interface.h:72
gvec_data gvec
datastruct holding global IMP state variable output
SF::vector< double > el_scale
optionally provided per-element params scale
Definition: fem_types.h:116
void dsp_from_cnt(const std::vector< T > &cnt, std::vector< T > &dsp)
Compute displacements from counts.
Definition: kdpart.hpp:140
void get_time(double &tm)
Definition: basics.h:436
Basic_physic * get_physics(physic_t p, bool error_if_missing)
Convinience function to get a physics.
Definition: sim_utils.cc:864
#define EXP_POSTPROCESS
Definition: sim_utils.h:161
int write_trace()
write traces to file
Definition: signals.h:701
T sum(const vector< T > &vec)
Compute sum of a vector&#39;s entries.
Definition: SF_vector.h:340
void open_trace(MULTI_IF *MIIF, int n_traceNodes, int *traceNodes, int *label, opencarp::sf_mesh *imesh)
Set up ionic model traces at some global node numbers.
int * N_Nodes
#nodes for each IMP
Definition: MULTI_ION_IF.h:200
sf_vec * old_vm
older Vm needed for 2nd order dT
Definition: electrics.h:101
void destroy()
Definition: ionics.cc:52
size_t root_count_ascii_lines(std::string file, MPI_Comm comm)
count the lines in a ascii file
int check_quiescence(double tm, double dt)
check for quiescence
Definition: electrics.cc:1641
V timing(V &t2, const V &t1)
Definition: basics.h:448
sig::time_trace wave
wave form of stimulus pulse
Definition: stimulate.h:96
void setup(int idx)
Setup from a param stimulus index.
Definition: stimulate.cc:161
void rebuild_matrices(MaterialType *mtype, limpet::MULTI_IF &miif, FILE_SPEC logger)
Definition: electrics.cc:1132
int regID
region ID
Definition: fem_types.h:95
File descriptor struct.
Definition: basics.h:133
void outer_prod(const Point &a, const Point &b, const double s, double *buff, const bool add=false)
Definition: SF_container.h:95
opencarp::sf_vec * gdata[NUM_IMP_DATA_TYPES]
data used by all IMPs
Definition: MULTI_ION_IF.h:216
void update_repolarization_times_from_rd(sf_vec &Vmv, sf_vec &Vmv_old, double time)
Checks if Vm crosses the repolarization threshold during each time step of the parabolic solver and u...
int set_dir(IO_t dest)
Definition: sim_utils.cc:483
void dup_IMP_node_state(IonIfBase &IF, int from, int to, GlobalData_t **localdata)
size_t size() const
The current size of the vector.
Definition: SF_vector.h:104
void log_stats(double tm, bool cflg)
Definition: timers.cc:27
int nsubregs
#subregions forming this region
Definition: fem_types.h:96
Electrical stimulation functions.
void log_msg(FILE_SPEC out, int level, unsigned char flag, const char *fmt,...)
Definition: basics.cc:72
#define NONL
Definition: basics.h:312
void init_sv_gvec(gvec_data &GVs, limpet::MULTI_IF *miif, sf_vec &tmpl, igb_output_manager &output_manager)
Definition: ionics.cc:566
int get_num_node() const
Gets the number of nodes handled by this IMP.
Definition: ION_IF.cc:148
bool mesh_is_registered(const mesh_t gt)
check wheter a SF mesh is set
Definition: sf_interface.cc:59
stim_t type
type of stimulus
Definition: stimulate.h:136
void save_eikonal_state(const char *tsav_ext)
Save the current state of variables related to the Eikonal simulation to a file to initialize a futur...
size_t l_numpts
local number of points
Definition: SF_container.h:389
lin_solver_stats stats
Definition: electrics.h:120
double x
Definition: SF_container.h:66
void output_initial_activations()
output one nodal vector of initial activation time
Definition: electrics.cc:1750
double y
Definition: SF_container.h:67
void dump_trace(MULTI_IF *MIIF, limpet::Real time)
size_t g_numpts
global number of points
Definition: SF_container.h:388
igb_output_manager output_manager_time
class handling the igb output
const T * begin() const
Pointer to the vector&#39;s start.
Definition: SF_vector.h:116
void f_close(FILE_SPEC &f)
Close a FILE_SPEC.
Definition: basics.cc:162
void compute_step()
Definition: ionics.cc:35
limpet::MULTI_IF * miif
Definition: ionics.h:66
sf_mat * mass_i
lumped for parabolic problem
Definition: electrics.h:109
sf_vec * IIon
ionic currents
Definition: electrics.h:98
void register_output(sf_vec *inp_data, const mesh_t inp_meshid, const int dpn, const char *name, const char *units, const SF::vector< mesh_int_t > *idx=NULL, bool elem_data=false)
Register a data vector for output.
Definition: sim_utils.cc:1447
void assign(InputIterator s, InputIterator e)
Assign a memory range.
Definition: SF_vector.h:161
const vector< T > & algebraic_nodes() const
Getter function for the local indices forming the local algebraic node set.
vector< S > she
sheet direction
Definition: SF_container.h:409
Diffusion Reaction Eikonal Alternant Model (DREAM) based on the electrics physics class...
size_t g_numelem
global number of elements
Definition: SF_container.h:386
double SF_real
Use the general double as real type.
Definition: SF_globals.h:38
Basic utility structs and functions, mostly IO related.
std::string timer_unit(const int timer_id)
figure out units of a signal linked to a given timer
void translate(int id)
convert legacy definitions to new format
Definition: stimulate.cc:100
stim_physics phys
physics of stimulus
Definition: stimulate.h:167
void apply_stim_to_vector(const stimulus &s, sf_vec &vec, bool add)
Definition: electrics.cc:428
SF::vector< RegionSpecs > regions
array with region params
Definition: fem_types.h:115
void init_logger(const char *filename)
void setup_dataout(const int dataout, std::string dataout_vtx, mesh_t grid, SF::vector< mesh_int_t > *&restr, bool async)
Definition: electrics.cc:588
void resize(size_t n)
Resize a vector.
Definition: SF_vector.h:209
vector< elem_t > type
element type
Definition: SF_container.h:406
int get_size(MPI_Comm comm=PETSC_COMM_WORLD)
Definition: basics.h:290
SF::scattering * register_scattering(const int from, const int to, const SF::SF_nbr nbr, const int dpn)
Register a scattering between to grids, or between algebraic and nodal representation of data on the ...
Definition: sf_interface.cc:65
void invert_3x3(S *ele, S &det)
Definition: dense_mat.hpp:452
int stimidx_from_timeridx(const SF::vector< stimulus > &stimuli, const int timer_id)
determine link between timer and stimulus
Definition: electrics.cc:747
void update_repolarization_times(const Ionics &ion)
Estimates initial repolarization times (T_R) in Step D of DREAM.
sf_mat * rhs_parab
rhs matrix to solve parabolic
Definition: electrics.h:110
centralize time managment and output triggering
Definition: timer_utils.h:73
int get_rank(MPI_Comm comm=PETSC_COMM_WORLD)
Definition: basics.h:276
T local_nodal_to_local_petsc(const meshdata< T, S > &mesh, int rank, T local_nodal)
void compute_bc()
Computes and updates boundary counditions in cycFIM.
grid_t
An electrics grid identifier to distinguish between intra and extra grids.
FILE_SPEC f_open(const char *fname, const char *mode)
Open a FILE_SPEC.
Definition: basics.cc:135
Container for a PETSc VecScatter.