ionicsOnFace.cc

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
// ----------------------------------------------------------------------------
 
#if WITH_EMI_MODEL
#include "ionicsOnFace.h"
 
#include "SF_init.h"
 
namespace opencarp {
 
void initialize_sv_dumps_onFace(limpet::MULTI_IF *pmiif, IMPregion_EMI* reg, int id, double t, double dump_dt);
 
namespace {
 
bool imp_region_emi_is_assigned(const IMPregion_EMI& region, int idx)
{
  if (!(region.im && strlen(region.im) > 0)) return false;
  if (idx < 2) return true;
  return region.num_IDs > 0;
}
 
int effective_num_imp_regions_emi()
{
  const int configured_nreg = param_globals::num_imp_regions;
  int nreg = configured_nreg;
  while (nreg > 2 && !imp_region_emi_is_assigned(param_globals::imp_region_emi[nreg - 1], nreg - 1)) {
    --nreg;
  }
 
  if (nreg < configured_nreg) {
    static bool warned_once = false;
    int rank = 0;
    MPI_Comm_rank(MPI_COMM_WORLD, &rank);
    if (!warned_once && rank == 0) {
      log_msg(NULL, 1, ECHO,
              "Warning: num_imp_regions=%d but only %d imp_region_emi entries are assigned; ignoring trailing unassigned regions.\n",
              configured_nreg, nreg);
      warned_once = true;
    }
  }
 
  return nreg;
}
 
}  // namespace
 
void IonicsOnFace::compute_step()
{
  double t1, t2;
  get_time(t1);
 
  miif->compute_ionic_current();
 
  double comp_time = timing(t2, t1);
  this->compute_time += comp_time;
 
  comp_stats.calls++;
  comp_stats.tot_time += comp_time;
 
  if(user_globals::tm_manager->trigger(iotm_console))
    comp_stats.log_stats(user_globals::tm_manager->time, false);
}
 
void IonicsOnFace::destroy()
{
  miif->free_MIIF();
}
 
void IonicsOnFace::output_step()
{}
 
void IonicsOnFace::initialize()
{
  double t1, t2;
  get_time(t1);
 
  set_dir(OUTPUT);
 
  // initialize generic logger for ODE timings per time_dt
  comp_stats.init_logger("ODE_stats.dat");
 
  double tstart = 0.;
  sf_mesh & mesh = get_mesh(ion_domain);
  int loc_size = mesh.l_numelem; 
 
  // create ionic current vector and register it
  sf_vec *IIon;
  sf_vec *Vmv;
  SF::init_vector(&IIon, mesh, 1, sf_vec::elemwise);
  SF::init_vector(&Vmv, mesh, 1, sf_vec::elemwise);
  register_data(IIon, iion_emi_itf);
  register_data(Vmv,  vm_emi_itf);
 
  // setup miif
  miif = new limpet::MULTI_IF();
  // store IIF_IDs and Plugins in arrays
  miif->name = "myocardium";
  miif->gdata[limpet::Vm]   = Vmv;
  miif->gdata[limpet::Iion] = IIon;
 
  const int num_imp_regions = effective_num_imp_regions_emi();
  SF::vector<RegionSpecs_EMI> rs(num_imp_regions);
 
  if (rs.size() <2) {
      log_msg(NULL, 5, ECHO, "\tError: num_imp_regions must be at least 2: the first region is reserved for the ionic model, and the second is designated for the gap junction.\n");
      log_msg(NULL, 5, ECHO, "set a default ionic model and gapjuntion model in parameter file\\n");
      exit(1);
    }
 
  // Region 0 and 1 are implicit defaults:
  // - region 0: default membrane faces
  // - region 1: default gap-junction faces
  rs[0].nsubregs = 0;
  rs[0].subregtags = nullptr;
  rs[1].nsubregs = 0;
  rs[1].subregtags = nullptr;
 
  // imp_region_emi[2+] assign models to specific face tag pairs like "t1:t2"
  for (size_t i = 2; i < rs.size(); i++ ) {
    rs[i].nsubregs   = param_globals::imp_region_emi[i].num_IDs;
    rs[i].subregtags = new std::string[rs[i].nsubregs]; // Allocate memory for the string array
    for (int j = 0; j < rs[i].nsubregs;j++) {
      std::string t1t2 = param_globals::imp_region_emi[i].ID[j];
      size_t colon_pos = t1t2.find(':');
      if (colon_pos == std::string::npos) {
        std::cerr << "Error: ':' the face tags are not defined properly in input file, it should be tag1:tag2 as a string" << std::endl;
        exit(1);
      }
 
      rs[i].subregtags[j] = param_globals::imp_region_emi[i].ID[j];  // Convert char* to std::string 
    }
  }
 
  SF::vector<int> reg_mask;
  region_mask_onFace(ion_domain, tags_data, line_face, tri_face, quad_face,
                     map_vertex_tag_to_dof, map_elem_uniqueFace_to_tags, intra_tags,
                     rs, reg_mask, true, "imp_region_emi");
 
  for (size_t i = 0; i < rs.size(); i++ ) {
    delete[] rs[i].subregtags;
    rs[i].subregtags = nullptr;
  }
 
  int purkfLen = 1;
  tstart = setup_MIIF(loc_size, num_imp_regions, param_globals::imp_region_emi,
                      reg_mask.data(), param_globals::start_statef, param_globals::num_adjustments,
                      param_globals::adjustment, param_globals::dt, purkfLen > 0);
 
  miif->extUpdateVm = !param_globals::operator_splitting;
 
  // if we start at a non-zero time (i.e. we have restored a state), we notify the
  // timer manager
  if(tstart > 0.) {
    log_msg(logger, 0, 0, "Changing simulation start time to %.2lf", tstart);
    user_globals::tm_manager->setup(param_globals::dt, tstart, param_globals::tend);
    user_globals::tm_manager->reset_timers();
  }
 
  set_dir(INPUT);
 
  this->initialize_time += timing(t2, t1);
}
 
double IonicsOnFace::setup_MIIF(int nnodes, int nreg, IMPregion_EMI* impreg, int* mask,
           const char *start_fn, int numadjust, IMPVariableAdjustment *adjust,
           double time_step, bool close)
{
  double tstart = 0;
 
  miif->N_IIF     = nreg;
  miif->numNode   = nnodes;
  miif->iontypes  = {};
  miif->numplugs  = (int*)calloc( miif->N_IIF, sizeof(int));
  miif->plugtypes = std::vector<limpet::IonTypeList>(miif->N_IIF);
  miif->targets = std::vector<limpet::Target>(miif->N_IIF, limpet::Target::AUTO);
 
  log_msg(logger,0,ECHO, "\nSetting up ionic models on EMI Face and plugins\n" \
                             "-----------------------------------\n\n" \
                             "Assigning IMPS to tagged regions:" );
 
  for (int i=0;i<miif->N_IIF;i++) {
    auto pT = limpet::get_ion_type(std::string(impreg[i].im));
    if (pT != NULL)
    {
      miif->iontypes.push_back(*pT);
      log_msg(logger, 0, ECHO|NONL, "\tIonic model: %s to tag region(s)", impreg[i].im);
 
      if(impreg[i].num_IDs > 0) {
        for(int j = 0; j < impreg[i].num_IDs; j++)
          log_msg(logger,0,ECHO|NONL, " [%s],", impreg[i].ID[j]);
        log_msg(logger,0,ECHO,"\b.");
      }
      else {
        log_msg(logger,0,ECHO, " [0] (implicitely)");
      }
    }
    else {
      log_msg(NULL,5,ECHO, "Illegal IM specified: %s\n", impreg[i].im );
      log_msg(NULL,5,ECHO, "Run bench --list-imps for a list of all available models.\n" );
      EXIT(1);
    }
    if (limpet::get_plug_flag( impreg[i].plugins, &miif->numplugs[i], miif->plugtypes[i]))
    {
      if(impreg[i].plugins[0] != '\0')  {
        log_msg(logger,0, ECHO|NONL, "\tPlug-in(s) : %s to tag region(s)",  impreg[i].plugins);
 
        for(int j = 0; j < impreg[i].num_IDs; j++)
          log_msg(logger,0,ECHO|NONL, " [%d],", impreg[i].ID[j]);
        log_msg(logger,0,ECHO,"\b.");
      }
    }
    else  {
      log_msg(NULL,5,ECHO,"Illegal plugin specified: %s\n", impreg[i].plugins);
      log_msg(NULL,5,ECHO, "Run bench --list-imps for a list of all available plugins.\n" );
      EXIT(1);
    }
  }
 
  miif->IIFmask = (limpet::IIF_Mask_t*)calloc(miif->numNode, sizeof(limpet::IIF_Mask_t));
 
  // The mask is a nodal vector with the region IDs of each node.
  // It is already reduced during the region_mask call to guarantee unique values
  // for overlapping interface nodes
  if (mask) {
    for (int i=0; i<miif->numNode; i++)
      miif->IIFmask[i] = (limpet::IIF_Mask_t) mask[i];
  }
 
  miif->initialize_MIIF();
 
  for (int i=0;i<miif->N_IIF;i++) {
    // the IMP tuning does not handle spaces well, thus we remove them here
    remove_char(impreg[i].im_param, strlen(impreg[i].im_param), ' ');
    miif->IIF[i]->tune(impreg[i].im_param, impreg[i].plugins, impreg[i].plug_param);
  }
 
  set_dir(INPUT);
  miif->initialize_currents(time_step, param_globals::ode_fac);
 
  // overriding initial values goes here
  // read in single cell state vector and spread it out over the entire region
  set_dir(INPUT);
  for (int i=0;i<miif->N_IIF;i++) {
    if (impreg[i].im_sv_init && strlen(impreg[i].im_sv_init) > 0)
      if (read_sv(miif, i, impreg[i].im_sv_init)) {
        log_msg(NULL, 5, ECHO|FLUSH, "State vector initialization failed for %s.\n", impreg[i].name);
        EXIT(-1);
      }
  }
 
  if( !start_fn || strlen(start_fn)>0 )
    tstart = (double) miif->restore_state(start_fn, ion_domain, close);
 
  for (int i=0; i<numadjust; i++)
  {
    set_dir(INPUT);
 
    SF::vector<int>    indices;
    SF::vector<double> values;
    bool restrict_to_algebraic = true;
 
    sf_mesh & imesh = get_mesh(ion_domain);
    std::map<std::string,std::string> metadata;
    read_metadata(adjust[i].file, metadata, PETSC_COMM_WORLD);
 
    SF::SF_nbr nbr = SF::NBR_REF;
    if(metadata.count("grid") && metadata["grid"].compare("intra") == 0) {
        nbr = SF::NBR_SUBMESH;
    }
    read_indices_with_data(indices, values, adjust[i].file, imesh, nbr, restrict_to_algebraic, 1, PETSC_COMM_WORLD);
 
    int rank = get_rank();
 
    for(size_t gi = 0; gi < indices.size(); gi++)
      indices[gi] = SF::local_nodal_to_local_petsc(imesh, rank, indices[gi]);
 
    // debug, output parameters on global intracellular vector
    if(adjust[i].dump) {
      sf_vec* adjPars;
      SF::init_vector(&adjPars, imesh, 1, sf_vec::algebraic);
      adjPars->set(indices, values);
 
      set_dir(OUTPUT);
      char fname[2085];
      snprintf(fname, sizeof fname, "adj_%s_perm.dat", adjust[i].variable);
      adjPars->write_ascii(fname, false);
 
      // get the scattering to the canonical permutation
      SF::scattering* sc = get_permutation(ion_domain, PETSC_TO_CANONICAL, 1);
      if(sc == NULL) {
        log_msg(0,3,0, "%s warning: PETSC_TO_CANONICAL permutation needed registering!", __func__);
        sc = register_permutation(ion_domain, PETSC_TO_CANONICAL, 1);
      }
 
      (*sc)(*adjPars, true);
      snprintf(fname, sizeof fname, "adj_%s_canonical.dat", adjust[i].variable);
      adjPars->write_ascii(fname, false);
    }
 
    int nc = miif->adjust_MIIF_variables(adjust[i].variable, indices, values);
    log_msg(logger, 0, 0, "Adjusted %d values for %s", nc, adjust[i].variable);
  }
 
  set_dir(OUTPUT);
 
  for (int i=0;i<miif->N_IIF;i++)
    initialize_sv_dumps_onFace(miif, impreg+i, i, tstart, param_globals::spacedt);
 
  return tstart;
}
 
void initialize_sv_dumps_onFace(limpet::MULTI_IF *pmiif, IMPregion_EMI* reg, int id, double t, double dump_dt)
{
  char  svs[1024], plgs[1024], plgsvs[1024], fname[1024];
 
  strcpy(svs, reg->im_sv_dumps ? reg->im_sv_dumps : "");
  strcpy(plgs, reg->plugins ? reg->plugins : "");
  strcpy(plgsvs, reg->plug_sv_dumps ? reg->plug_sv_dumps : "");
 
  if( !(strlen(svs)+strlen(plgsvs) ) )
    return;
 
  /* The string passed to the "reg_name" argument (#4) of the sv_dump_add
  *     function is supposed to be "region name". It's only purpose is to
  *     provide the base name for the SV dump file, eg: Purkinje.Ca_i.bin.
  *     Thus, we pass:  [vofile].[reg name], Otherwise, dumping SVs in
  *     batched runs would be extremely tedious.
  */
  strcpy(fname, param_globals::vofile);                              // [vofile].igb.gz
 
  if( !reg->name ) {
    log_msg(NULL, 5, ECHO, "%s: a region name must be specified\n", __func__ );
    exit(0);
  }
 
  // We want to convert vofile.igb or vofile.igb.gz to vofile.regname
  char* ext_start = NULL;
  ext_start = strstr(fname, "igb.gz");
  if(ext_start == NULL) ext_start = strstr(fname, "igb");
  if(ext_start == NULL) ext_start = fname + strlen(fname);
  strcpy(ext_start, reg->name);
 
  pmiif->sv_dump_add_by_name_list(id, reg->im, fname, svs, plgs, plgsvs, t, dump_dt);
}
 
bool check_tags_in_elems_onFace(std::vector<std::string> & tags_data, SF::vector<RegionSpecs_EMI> & regspec,
                         const char* gridname, const char* reglist)
{
  bool AllTagsExist = true;
  hashmap::unordered_set<std::string> tagset;
 
  tagset.insert(tags_data.begin(), tags_data.end());
 
  // cycle through all user-specified regions
  for (size_t reg=0; reg<regspec.size(); reg++)
    // cycle through all tags which belong to the region
    for (int k=0; k<regspec[reg].nsubregs; k++) {
      // check whether this tag exists in element list
      int n = 0;
      size_t colon_pos = regspec[reg].subregtags[k].find(':');
      int tag1 = std::stoi(regspec[reg].subregtags[k].substr(0, colon_pos));
      int tag2 = std::stoi(regspec[reg].subregtags[k].substr(colon_pos + 1));
 
      // considered the other combinations of tags
      std::string inverse_pair;
      inverse_pair = std::to_string(tag2) + ":" + std::to_string(tag1);
      if(tagset.count(regspec[reg].subregtags[k]) || tagset.count(inverse_pair)) {
        n++;
      }
 
      // globalize n
      int N = get_global(n, MPI_SUM);
 
      if (N==0) {
        log_msg(NULL, 3, ECHO,
                "on face Region tag %s in %s[%d] not found in element list for %s grid.\n",
                regspec[reg].subregtags[k].c_str(), reglist, reg, gridname);
        AllTagsExist = false;
      }
    }
 
  if (!AllTagsExist) {
    log_msg(NULL, 4, ECHO,"Assigned wrong pair of tags on the face of EMI surface mesh!\n"
                          "Check region ID specs in input files!\n");
  }
 
  return AllTagsExist;
}
 
inline bool pair_is_gapjunction(const std::pair<mesh_int_t, mesh_int_t>& tags,
                                const hashmap::unordered_set<int>& intra_tags)
{
  return intra_tags.find(tags.first) != intra_tags.end() &&
         intra_tags.find(tags.second) != intra_tags.end();
}
 
void region_mask_onFace(mesh_t meshspec, 
                        std::vector<std::string> & tags_data,                   
                        hashmap::unordered_map<SF::tuple<mesh_int_t>, 
                        std::pair<SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>,
                        SF::emi_face<mesh_int_t,SF::tuple<mesh_int_t>>>> & line_face,
                        hashmap::unordered_map<SF::triple<mesh_int_t>, 
                        std::pair<SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>,
                        SF::emi_face<mesh_int_t,SF::triple<mesh_int_t>>>> & tri_face,
                        hashmap::unordered_map<SF::quadruple<mesh_int_t>, 
                        std::pair<SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>, 
                        SF::emi_face<mesh_int_t,SF::quadruple<mesh_int_t>>>> & quad_face,
                        hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, mesh_int_t> & map_vertex_tag_to_dof,
                        hashmap::unordered_map<mesh_int_t, std::pair<mesh_int_t, mesh_int_t>> & map_elem_uniqueFace_to_tags,
                        const hashmap::unordered_set<int> & intra_tags,
                        SF::vector<RegionSpecs_EMI> & regspec,
                        SF::vector<int> & regionIDs, bool mask_elem, const char* reglist)
{
  if(regspec.size() == 1) return;
 
  sf_mesh & mesh = get_mesh(meshspec);
  SF::element_view<mesh_int_t,mesh_real_t> eview(mesh, SF::NBR_PETSC); // error is here!!!
 
  // initialize the list with the default regionID, 0
  size_t rIDsize = mask_elem ? mesh.l_numelem : mesh.l_numpts;
 
  size_t nelem = mesh.l_numelem;
 
  // check whether all specified tags exist in the element list
  check_tags_in_elems_onFace(tags_data, regspec, mesh.name.c_str(), reglist);
 
  regionIDs.assign(rIDsize, 0);
 
  int* rid = regionIDs.data();
  hashmap::unordered_map<std::string,int> tag_to_reg;
 
  // we generate a map from tags to region IDs. This has many benefits, mainly we can check
  // whether a tag is assigned to multiple regions and simplify our regionIDs filling loop
  for (size_t reg=0; reg < regspec.size(); reg++) {
    int err = 0;
    for (int k=0; k < regspec[reg].nsubregs; k++) {
      std::string curtag = regspec[reg].subregtags[k];
      if(tag_to_reg.count(curtag)) err++;
 
      if(get_global(err, MPI_SUM))
        log_msg(0,4,0, "%s warning: Tag idx %s is assigned to multiple regions!\n"
                       "Its final assignment will be to the highest assigned region ID!",
                __func__, curtag.c_str());
 
      tag_to_reg[curtag] = reg;
    }
  }
 
  std::vector<int> mx_tag( rIDsize, -1 );
 
  // cycle through the element list
  for(size_t eidx=0; eidx<nelem; eidx++)
  {
    std::vector<int> elem_nodes;
    mesh_int_t tag = mesh.tag[eidx];
    std::string result_pair_orginal;
    std::string result_pair_reverse;
    std::pair<mesh_int_t,mesh_int_t> value = map_elem_uniqueFace_to_tags[eidx];
    result_pair_orginal = std::to_string(value.first) + ":" + std::to_string(value.second);
    result_pair_reverse = std::to_string(value.second) + ":" + std::to_string(value.first);
    rid[eidx] = pair_is_gapjunction(value, intra_tags) ? 1 : 0;
 
    if (tag_to_reg.count(result_pair_orginal)) {
      rid[eidx] = tag_to_reg[result_pair_orginal];
    } else if(tag_to_reg.count(result_pair_reverse)){
      rid[eidx] = tag_to_reg[result_pair_reverse];
    }
    if(tag!=value.first){
      log_msg(NULL, 5, ECHO, "error in tags on surface mesh with the unique faces!!!.\\n");
      exit(1);
    }
  }
}
 
double IonicsOnFace::timer_val(const int timer_id)
{
  double val = std::nan("NaN");
  return val;
}
 
std::string IonicsOnFace::timer_unit(const int timer_id)
{
  std::string s_unit;
  return s_unit;
}
 
void compute_IIF_OnFace(limpet::IonIfBase& pIF, limpet::GlobalData_t** impdata, int n)
{
  if (impdata[limpet::Iion] != NULL)
    impdata[limpet::Iion][n] = 0;
 
  pIF.for_each([&](limpet::IonIfBase& imp) {
    update_ts(&imp.get_tstp());
    imp.compute(n, n + 1, impdata);
  });
}
 
void* find_SV_in_IMP_onFace(limpet::MULTI_IF* miif, const int idx, const char *IMP, const char *SV,
                     int* offset, int* sz)
{
  if(strcmp(IMP, miif->iontypes[idx].get().get_name().c_str()) == 0) {
    return (void*) miif->iontypes[idx].get().get_sv_offset(SV, offset, sz);
  }
  else {
    for( int k=0; k<miif->numplugs[idx]; k++ )
      if(strcmp(IMP, miif->plugtypes[idx][k].get().get_name().c_str()) == 0) {
        return (void*) miif->plugtypes[idx][k].get().get_sv_offset(SV, offset, sz);
      }
  }
 
  return NULL;
}
 
 
void alloc_gvec_data_onFace(const int nGVcs, const int nRegs, const int nEmiRegs,
                     GVecs *prmGVecs, gvec_data_OnFace &glob_vecs)
{
  glob_vecs.nRegs = nRegs;
 
  if (nGVcs) {
    glob_vecs.vecs.resize(nGVcs);
 
    for (size_t i = 0; i < glob_vecs.vecs.size(); i++) {
      sv_data_onFace &gvec = glob_vecs.vecs[i];
 
      gvec.name    = dupstr(prmGVecs[i].name);
      gvec.units   = dupstr(prmGVecs[i].units);
      gvec.bogus   = prmGVecs[i].bogus;
 
      gvec.imps    = (char**) calloc(nRegs, sizeof(char *));
      gvec.svNames = (char**) calloc(nRegs, sizeof(char *));
      gvec.svSizes = (int*)   calloc(nRegs, sizeof(int));
      gvec.svOff   = (int*)   calloc(nRegs, sizeof(int));
      gvec.getsv   = (void**) calloc(nRegs, sizeof(limpet::SVgetfcn));
 
      for (int j = 0; j < nRegs; j++) {
        if (strlen(prmGVecs[i].imp)) gvec.imps[j] = dupstr(prmGVecs[i].imp);
#ifdef WITH_PURK
        else if (j < nEmiRegs)
          gvec.imps[j] = dupstr(param_globals::imp_region_emi[j].im);
        else
          gvec.imps[j] = dupstr(param_globals::PurkIon[j - nEmiRegs].im);
#else
        else if (j < nEmiRegs)
          gvec.imps[j] = dupstr(param_globals::imp_region_emi[j].im);
#endif
        gvec.svNames[j] = dupstr(prmGVecs[i].ID[j]);
      }
    }
  }
}
 
void init_sv_gvec_onFace(gvec_data_OnFace& GVs, limpet::MULTI_IF* miif, sf_vec & tmpl,
                  igb_output_manager & output_manager)
{
  GVs.inclPS = false;
  // int num_purk_regions = GVs->inclPS ? purk->ion.N_IIF : 0;
  int num_purk_regions = 0;
  int nEmiRegs         = effective_num_imp_regions_emi();
  int nRegs            = nEmiRegs + num_purk_regions;
 
  alloc_gvec_data_onFace(param_globals::num_gvecs, nRegs, nEmiRegs, param_globals::gvec, GVs);
 
#ifdef WITH_PURK
  if (GVs->inclPS) sample_PS_ionSVs(purk);
#endif
 
  for (unsigned int i = 0; i < GVs.vecs.size(); i++) {
    sv_data_onFace & gv = GVs.vecs[i];
    int noSV    = 0;
 
    for (int j = 0; j < miif->N_IIF; j++) {
      gv.getsv[j] = find_SV_in_IMP_onFace(miif, j, gv.imps[j], gv.svNames[j],
                                   gv.svOff + j, gv.svSizes + j);
 
      if (gv.getsv[j] == NULL) {
        log_msg(NULL, 3, ECHO, "\tWarning: SV(%s) not found in region %d\n", gv.svNames[j], j);
        noSV++;
      }
    }
 
    SF::init_vector(&gv.ordered, &tmpl);
    output_manager.register_output(gv.ordered, intra_elec_msh, 1, gv.name, gv.units);
 
#ifdef WITH_PURK
    // same procedure for Purkinje
    if (GVs->inclPS) {
      IF_PURK_PROC(purk) {
        MULTI_IF* pmiif = &purk->ion;
        for (int j = miif->N_IIF; j < nRegs; j++) {
          gv.getsv[j] = find_SV_in_IMP_onFace(pmiif, j - miif->N_IIF, gv.imps[j],
                                       gv.svNames[j], gv.svOff + j, gv.svSizes + j);
 
          if (gv.getsv[j] == NULL) {
            LOG_MSG(NULL, 3, ECHO, "\tWarning: state variable \"%s\" not found in region %d\n", gv.svNames[j], j);
            noSV++;
          }
        }
        RVector_dup(purk->vm_pt, &gv.orderedPS_PS);
        MYO_COMM(purk);
      }
      RVector_dup(purk->vm_pt_over, &gv.orderedPS);
      initialize_grid_output(grid, NULL, tmo, intra_elec_msh, GRID_WRITE, 0., 1., gv.units, gv.orderedPS,
                             1, gv.GVcName, -purk->npt, param_globals::output_level);
    }
#endif
 
    if (noSV == nRegs) {
      log_msg(NULL, 5, ECHO, "\tError: no state variables found for global vector %d\n", i);
      log_msg(NULL, 5, ECHO, "Run bench --imp=YourModel --imp-info to get a list of all parameters.\\n");
      exit(1);
    }
  }
}
 
void assemble_sv_gvec_onFace(gvec_data_OnFace & gvecs, limpet::MULTI_IF *miif)
{
  for(size_t i=0; i<gvecs.vecs.size(); i++ ) {
    sv_data_onFace & gv = gvecs.vecs[i];
 
    // set to the defalt value
    gv.ordered->set(gv.bogus);
    SF_int start, stop;
    gv.ordered->get_ownership_range(start, stop);
 
    for( int n = 0; n<miif->N_IIF; n++ ) {
      if( !gv.getsv[n] ) continue;
 
      SF::vector<SF_real> data   (miif->N_Nodes[n]);
      SF::vector<SF_int>  indices(miif->N_Nodes[n]);
 
      for( int j=0; j<miif->N_Nodes[n]; j++ ) {
        indices[j] = miif->NodeLists[n][j] + start;
        data[j]    = ((limpet::SVgetfcn)(gv.getsv[n]))( *miif->IIF[n], j, gv.svOff[n]);
      }
 
      bool add = false;
      gv.ordered->set(indices, data, add);
    }
 
#ifdef WITH_PURK
    // include purkinje in sv dump, if exists
    if(gvecs->inclPS) {
      RVector_set( gv.orderedPS, gv.bogus ); // set this to the default value
 
      if(IS_PURK_PROC(purk))
      {
        MULTI_IF *pmiif = &purk->ion;
        Real *data = new Real[purk->gvec_cab.nitems];
 
        int ci=0;
        for( int n = 0; n<pmiif->N_IIF; n++ ) {
          int gvidx = n+miif->N_IIF;
          if( !gv.getsv[gvidx] ) continue;
          for( int j=0; j<purk->gvec_ion[n].nitems; j++ )
            data[ci++] = ((SVgetfcn)(gv.getsv[gvidx]))( pmiif->IIF+n, ((int*)(purk->gvec_ion[n].data))[j],
                                                                                          gv.svOff[gvidx]);
        }
        RVector_setvals(gv.orderedPS, purk->gvec_cab.nitems, (int*)purk->gvec_cab.data, data, true);
        delete [] data;
      }
      RVector_sync( gv.orderedPS );
    }
#endif
  }
}
 
 
} // namespace opencarp
#endif