SF_fem_utils_emi.h

Go to the documentation of this file.

// ----------------------------------------------------------------------------
// openCARP is an open cardiac electrophysiology simulator.
//
// Copyright (C) 2020 openCARP project
//
// This program is licensed under the openCARP Academic Public License (APL)
// v1.0: You can use and redistribute it and/or modify it in non-commercial
// academic environments under the terms of APL as published by the openCARP
// project v1.0, or (at your option) any later version. Commercial use requires
// a commercial license (info@opencarp.org).
//
// This program is distributed without any warranty; see the openCARP APL for
// more details.
//
// You should have received a copy of the openCARP APL along with this program
// and can find it online: http://www.opencarp.org/license
// ----------------------------------------------------------------------------
 
#if WITH_EMI_MODEL
 
#ifndef _SF_FEM_EMI_H
#define _SF_FEM_EMI_H
 
#include "SF_container.h"
#include "SF_fem_utils.h"
namespace SF {
 
inline double computeTriangleArea(const SF::Point& p1, const SF::Point& p2, const SF::Point& p3) 
{
  // Define two edges of the triangle
  Point edge1 = p2 - p1;
  Point edge2 = p3 - p1;
 
  // Compute the cross product
  Point crossProduct =  cross(edge1,edge2);
 
  // Compute the magnitude of the cross product
  double crossProductMagnitude = mag(crossProduct);
 
  return 0.5 * crossProductMagnitude;
}
 
inline void compute_barycentric_coordinates_coefficients(const SF::Point& p1, const SF::Point& p2, const SF::Point& P, vector<double>& interpolationCoefficients, double & lengthSquared) 
{
  SF::Point d = p2 - p1;
  lengthSquared = inner_prod(d, d);
 
  SF::Point v = P - p1;
  double lambda2 = inner_prod(v, d) / lengthSquared;
  double lambda1 = 1.0 - lambda2;
 
  interpolationCoefficients[0] = lambda1;
  interpolationCoefficients[1] = lambda2;
}
 
inline void compute_barycentric_coordinates_coefficients(const SF::Point& p1, const SF::Point& p2, const SF::Point& p3, const SF::Point& P, vector<double>& interpolationCoefficients, double & area) 
{
    
    area = computeTriangleArea(p1, p2, p3);
    
    // Vectors for the edges of the triangle
    SF::Point v0 = p2 - p1;
    SF::Point v1 = p3 - p1;
    SF::Point v2 = P - p1;
 
    // Dot products
    double d00 = inner_prod(v0,v0);
    double d01 = inner_prod(v0,v1);
    double d11 = inner_prod(v1,v1);
    double d20 = inner_prod(v2,v0);
    double d21 = inner_prod(v2,v1);
 
    // Compute barycentric coordinates
    double denom = d00 * d11 - d01 * d01;
    double beta = (d11 * d20 - d01 * d21) / denom;
    double gamma = (d00 * d21 - d01 * d20) / denom;
    double alpha = 1.0 - beta - gamma;
    interpolationCoefficients[0] = alpha; // α
    interpolationCoefficients[1] = beta;  // β
    interpolationCoefficients[2] = gamma; // ζ
}
 
 
inline void compute_barycentric_coordinates_coefficients(const SF::Point& p1, const SF::Point& p2, const SF::Point& p3, const SF::Point& p4, const SF::Point& P, vector<double>& interpolationCoefficients, double & area) 
{    
  // Area: compute as sum of two triangles
  double area1 = computeTriangleArea(p1, p2, p4);
  double area2 = computeTriangleArea(p2, p3, p4);
  area = area1 + area2;
 
  // Define local axes
  SF::Point origin = p1;
  SF::Point u_dir = p2 - p1;
  SF::Point v_dir = p4 - p1;
 
  // Project P into the local (u,v) frame
  SF::Point d = P - origin;
  double u = inner_prod(d, u_dir) / inner_prod(u_dir, u_dir);
  double v = inner_prod(d, v_dir) / inner_prod(v_dir, v_dir);
 
  // Clamp u and v to [0,1] to avoid extrapolation (optional)
  u = SF::clamp(u, 0.0, 1.0);
  v = SF::clamp(v, 0.0, 1.0);
 
  // Bilinear shape functions
  interpolationCoefficients[0] = (1 - u) * (1 - v); // α
  interpolationCoefficients[1] = u * (1 - v);       // β
  interpolationCoefficients[2] = u * v;             // γ
  interpolationCoefficients[3] = (1 - u) * v;       // ζ
}
 
template<class S>
inline void compute_integrate_matrix_barycentric(vector<SF::Point> face_coordinates, SF_int nnodes, dmat<SF_real> & ebuff, dmat<SF_real> & ebuff_s, dmat<SF_real> & ebuff_counter, S mass_scale)
{
  //compute the barycentric coordinates: calculateCentroid
  double x = 0;
  double y = 0;
  double z = 0;
  for (int i = 0; i < nnodes; ++i)
  {
    x +=  face_coordinates[i].x;
    y +=  face_coordinates[i].y;
    z +=  face_coordinates[i].z;
  }
  SF::Point b; b.x = x/nnodes; b.y = y/nnodes; b.z = z/nnodes; 
 
  vector<double> interpolationCoefficients(nnodes);
  double area, weighted_area;
  if(nnodes==2){
    compute_barycentric_coordinates_coefficients(face_coordinates[0], face_coordinates[1], b, interpolationCoefficients, area);
  }
  else if(nnodes==3)
  {
    compute_barycentric_coordinates_coefficients(face_coordinates[0], face_coordinates[1], face_coordinates[2], b, interpolationCoefficients, area);
  }
  else if(nnodes==4){
    compute_barycentric_coordinates_coefficients(face_coordinates[0], face_coordinates[1], face_coordinates[2], face_coordinates[3], b, interpolationCoefficients, area);
  }
 
  weighted_area = area/nnodes;
 
  if(nnodes==2){
    ebuff[0][0] = interpolationCoefficients[0];
    ebuff[0][1] = interpolationCoefficients[1];
 
    ebuff_counter[0][0] = interpolationCoefficients[0];
    ebuff_counter[0][1] = interpolationCoefficients[1];
 
    ebuff_s[0][0] = mass_scale*weighted_area;
    ebuff_s[0][1] = mass_scale*weighted_area;
  }
  else if(nnodes==3){
    ebuff[0][0] = interpolationCoefficients[0];
    ebuff[0][1] = interpolationCoefficients[1];
    ebuff[0][2] = interpolationCoefficients[2];
 
    ebuff_counter[0][0] = interpolationCoefficients[0];
    ebuff_counter[0][1] = interpolationCoefficients[1];
    ebuff_counter[0][2] = interpolationCoefficients[2];
 
    ebuff_s[0][0] = mass_scale*weighted_area;
    ebuff_s[0][1] = mass_scale*weighted_area;
    ebuff_s[0][2] = mass_scale*weighted_area;
  }
  else if(nnodes==4){
    ebuff[0][0] = interpolationCoefficients[0];
    ebuff[0][1] = interpolationCoefficients[1];
    ebuff[0][2] = interpolationCoefficients[2];
    ebuff[0][3] = interpolationCoefficients[3];
 
    ebuff_counter[0][0] = interpolationCoefficients[0];
    ebuff_counter[0][1] = interpolationCoefficients[1];
    ebuff_counter[0][2] = interpolationCoefficients[2];
    ebuff_counter[0][3] = interpolationCoefficients[3];
 
    ebuff_s[0][0] = mass_scale*weighted_area;
    ebuff_s[0][1] = mass_scale*weighted_area;
    ebuff_s[0][2] = mass_scale*weighted_area;
    ebuff_s[0][3] = mass_scale*weighted_area;
  }
}
 
template<class tuple_key, class tuple_value, class tri_key, class tri_value, class quad_key, class quad_value, class T, class S>
inline void assemble_restrict_operator( abstract_matrix<T,S> & B,
                                        abstract_matrix<T,S> & Bi,
                                        abstract_matrix<T,S> & BsM,
                                        SF::vector<mesh_int_t> elemTag_surface_mesh, 
                                        hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, 
                                        std::pair<mesh_int_t,mesh_int_t>> map_vertex_tag_to_dof_petsc,
                                        hashmap::unordered_map<tuple_key,
                                        std::pair<tuple_value, tuple_value>> line_face,
                                        hashmap::unordered_map<tri_key,
                                        std::pair<tri_value, tri_value>> tri_face, 
                                        hashmap::unordered_map<quad_key,
                                        std::pair<quad_value, quad_value>> quad_face,
                                        const meshdata<mesh_int_t,mesh_real_t> & surface_mesh,
                                        const meshdata<mesh_int_t,mesh_real_t> & emi_mesh,
                                        S mass_scale)
{
  // we want to make sure that the element integrator fits the matrix
  T row_dpn = 1; T col_dpn = 1;
 
  const vector<SF_int> & petsc_nbr = emi_mesh.get_numbering(NBR_PETSC);// remove it
  const SF::vector<mesh_int_t> & rnod_emi = emi_mesh.get_numbering(SF::NBR_REF);
  hashmap::unordered_map<T,T> g2l_emi;
  hashmap::unordered_map<T,T> l2g_emi;
  for(size_t i=0; i<rnod_emi.size(); i++){
    g2l_emi[rnod_emi[i]] = i;
    l2g_emi[i] = rnod_emi[i];
  }
 
 
  const SF::vector<mesh_int_t> & rnod = surface_mesh.get_numbering(SF::NBR_REF);
  hashmap::unordered_map<T,T> g2l;
  hashmap::unordered_map<T,T> l2g;
 
  for(size_t i=0; i<rnod.size(); i++){
    g2l[rnod[i]] = i;
    l2g[i] = rnod[i];
  }
 
  // allocate row / col index buffers
  vector<SF_int> row_idx(SF_MAX_ELEM_NODES * row_dpn), row_idx_counter(SF_MAX_ELEM_NODES * row_dpn), col_idx(SF_MAX_ELEM_NODES * col_dpn), col_idx_counter(SF_MAX_ELEM_NODES * col_dpn);
 
  const SF::vector<mesh_int_t> & emi_surfmesh_elem =  surface_mesh.get_numbering(SF::NBR_ELEM_REF);
 
  // allocate elem index buffer
  dmat<SF_real> ebuff(SF_MAX_ELEM_NODES * row_dpn, SF_MAX_ELEM_NODES * col_dpn);
  dmat<SF_real> ebuff_s(SF_MAX_ELEM_NODES * row_dpn, SF_MAX_ELEM_NODES * col_dpn);
  dmat<SF_real> ebuff_counter(SF_MAX_ELEM_NODES * row_dpn, SF_MAX_ELEM_NODES * col_dpn);
 
  const T* con_emi = emi_mesh.con.data();
  const T* con = surface_mesh.con.data();
 
  // start with assembly
  for(size_t eidx=0; eidx < surface_mesh.l_numelem; eidx++)
  {
    T tag = surface_mesh.tag[eidx];
 
    // get the counter part
    std::vector<int> elem_nodes;
    std::vector<int> elem_nodes_old;  
    std::vector<int> petsc_first;
    std::vector<int> petsc_second;
 
    T tag_first = 0;
    T tag_second = 0;
 
    T rank_first = 0;
    T rank_second = 0;
 
    T mem_first = 0;
    T mem_second = 0;
 
    T sign = 1.0;
 
    for (int n = surface_mesh.dsp[eidx]; n < surface_mesh.dsp[eidx+1];n++)
    {
      T l_idx = surface_mesh.con[n];
 
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(l2g[l_idx],tag);
      mesh_int_t Index_new = map_vertex_tag_to_dof_petsc[Index_tag_old].first;
      elem_nodes.push_back(Index_new);
      elem_nodes_old.push_back(l2g[l_idx]);
    }
 
    std::sort(elem_nodes.begin(),elem_nodes.end());
    if(elem_nodes.size()==2){
      tuple_key key;
 
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      std::pair<tuple_value, tuple_value> value = line_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
 
      rank_first = value.first.rank;
      rank_second = value.second.rank; 
 
      mem_first = value.first.mem;
      mem_second = value.second.mem;
    }
    else if(elem_nodes.size()==3){
      tri_key key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      std::pair<tri_value, tri_value> value = tri_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
 
      rank_first = value.first.rank;
      rank_second = value.second.rank;
 
      mem_first = value.first.mem;
      mem_second = value.second.mem;
    }
    else if(elem_nodes.size()==4){
      quad_key key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      key.v4 = elem_nodes[3];
      std::pair<quad_value, quad_value> value = quad_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
 
      rank_first = value.first.rank;
      rank_second = value.second.rank;
 
      mem_first = value.first.mem;
      mem_second = value.second.mem;
    }
 
    if(tag != tag_first)  std::swap(tag_first, tag_second); 
 
    std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
 
    for (int indx = 0; indx < elem_nodes_old.size(); ++indx)
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old_first;
      Index_tag_old_first = std::make_pair(elem_nodes_old[indx],tag_first);
      
      auto it_first = map_vertex_tag_to_dof_petsc.find(Index_tag_old_first);
      if (it_first == map_vertex_tag_to_dof_petsc.end()) {
        std::cerr << "ERROR: tag_first=" << tag_first << " not found for vertex=" << elem_nodes_old[indx] << std::endl;
        throw std::runtime_error("PETSc index not found for first tag");
      }
      std::pair <mesh_int_t,mesh_int_t> newIndex_petsc_first = it_first->second;
      mesh_int_t newIndex_first = newIndex_petsc_first.first;
      petsc_first.push_back(newIndex_petsc_first.second);
    
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old_second;
      Index_tag_old_second = std::make_pair(elem_nodes_old[indx],tag_second);
      
      auto it_second = map_vertex_tag_to_dof_petsc.find(Index_tag_old_second);
      if (it_second == map_vertex_tag_to_dof_petsc.end()) {
        std::cerr << "ERROR: tag_second=" << tag_second << " not found for vertex=" << elem_nodes_old[indx] << std::endl;
        throw std::runtime_error("PETSc index not found for counter-face tag");
      }
      std::pair <mesh_int_t,mesh_int_t> newIndex_petsc_second = it_second->second;
      mesh_int_t newIndex_second = newIndex_petsc_second.first;
      petsc_second.push_back(newIndex_petsc_second.second);
 
      // elemTag_surface_mesh[eidx] == 1 means that the tag belongs to the extracellular region
      if(mem_first==1 and elemTag_surface_mesh[eidx]==1)   // condition for the membrane which the sign is always negative on extracellular side
        sign = -1.0;
      else if (mem_first==2 and tag<tag_second)    // condition for the gap juntion which the sign is always negative for smaller tag
        sign = -1.0;
    }
 
    SF_int nnodes = elem_nodes.size();
    vector<SF::Point> face_coordinates(nnodes);
    // calculate row/col indices of entries
    row_idx.resize(nnodes*row_dpn);
    row_idx_counter.resize(nnodes*row_dpn);
    col_idx.resize(nnodes*col_dpn);
    col_idx_counter.resize(nnodes*col_dpn);
 
    {
      for(SF_int i=0; i<nnodes; i++){
        // all the row sets with element index 
        for(short j=0; j<row_dpn; j++){
          row_idx.data()[i*row_dpn + j] = emi_surfmesh_elem[eidx]*row_dpn; 
        }
 
        for(short j=0; j<col_dpn; j++){
          col_idx[i*col_dpn + j] = (petsc_first[i])*col_dpn + j;
          col_idx_counter[i*col_dpn + j] = (petsc_second[i])*col_dpn + j;  
        }
      }
    }
 
    // assign ebuff and ebuff_s
    ebuff.assign(nnodes, nnodes, 0.0);
    ebuff_s.assign(nnodes, nnodes, 0.0);
 
    ebuff_counter.assign(nnodes, nnodes, 0.0);
 
    for (int n = surface_mesh.dsp[eidx], i = 0; n < surface_mesh.dsp[eidx+1];n++,i++)
    {
      T l_idx = surface_mesh.con[n];
 
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(l2g[l_idx],tag);
      mesh_int_t Index_new = map_vertex_tag_to_dof_petsc[Index_tag_old].first;
 
      double x = emi_mesh.xyz[g2l_emi[Index_new]*3+0];
      double y = emi_mesh.xyz[g2l_emi[Index_new]*3+1];
      double z = emi_mesh.xyz[g2l_emi[Index_new]*3+2];
 
      face_coordinates[i].x = x;
      face_coordinates[i].y = y;
      face_coordinates[i].z = z;
    }
    compute_integrate_matrix_barycentric(face_coordinates, nnodes, ebuff, ebuff_s, ebuff_counter, mass_scale);
 
    // add values into system matrix
    const bool add = true;
 
    // we simply build Bi in the same way as B, but without the sign change, so we end up with Bi*Iij_stim = Ib_stim = Ii + Ij
    Bi.set_values(row_idx, col_idx, ebuff.data(), add);
    Bi.set_values(row_idx, col_idx_counter, ebuff_counter.data(), add);
    // the operator matrix in order to project potential of dofs to the barycenteric of each face as transmembrane (vm = ui - uj)
    ebuff*=(sign);
    B.set_values(row_idx, col_idx, ebuff.data(), add);
    ebuff_counter*=(-sign);
    B.set_values(row_idx, col_idx_counter, ebuff_counter.data(), add); // counter face
 
    // the operator matrix in order to project the current density to the dofs of the face
    ebuff_s*=(sign);
    BsM.set_values(col_idx, row_idx, ebuff_s.data(), add);
  }
  // finish assembly and progress output
  B.finish_assembly();
  Bi.finish_assembly();
  BsM.finish_assembly();
}
 
template<class tuple_key, class tuple_value, class tri_key, class tri_value, class quad_key, class quad_value, class T, class S>
inline void assemble_lhs_emi(abstract_matrix<T,S> & mat,
                             abstract_matrix<T,S> & mat_surf,
                             const meshdata<mesh_int_t,mesh_real_t> & emi_mesh,
                             const meshdata<mesh_int_t,mesh_real_t> & surface_mesh,
                             hashmap::unordered_map<std::pair<mesh_int_t,mesh_int_t>, std::pair<mesh_int_t,mesh_int_t>> map_vertex_tag_to_dof_petsc,
                             hashmap::unordered_map<tuple_key,
                             std::pair<tuple_value, tuple_value>> line_face,
                             hashmap::unordered_map<tri_key,
                             std::pair<tri_value, tri_value>> tri_face, 
                             hashmap::unordered_map<quad_key,
                             std::pair<quad_value, quad_value>> quad_face,
                             matrix_integrator<mesh_int_t,mesh_real_t> & stiffness_integrator,
                             matrix_integrator<mesh_int_t, mesh_real_t> & mass_integrator,
                             S stiffness_scale,
                             S mass_scale)
{
  // we want to make sure that the element integrator fits the matrix
  T row_dpn = 1; T col_dpn = 1;
 
  // allocate row / col index buffers
  vector<SF_int> idx(SF_MAX_ELEM_NODES * col_dpn), idx_counter(SF_MAX_ELEM_NODES * col_dpn);
  vector<T> row_idx(SF_MAX_ELEM_NODES * row_dpn), col_idx(SF_MAX_ELEM_NODES * col_dpn);
 
  // allocate elem index buffer
  dmat<SF_real> ebuff(SF_MAX_ELEM_NODES * col_dpn, SF_MAX_ELEM_NODES * col_dpn);
  dmat<SF_real> ebuff_mass(SF_MAX_ELEM_NODES * col_dpn, SF_MAX_ELEM_NODES * col_dpn);
  // Assemble stiffness matrix (K)
  const vector<mesh_int_t> & stiffness_petsc_nbr = emi_mesh.get_numbering(NBR_PETSC);
  element_view<mesh_int_t, mesh_real_t> stiffness_view(emi_mesh, NBR_PETSC);
 
  // Assembly of global stiffness matrix into lhs matrix 
  for(size_t eidx=0; eidx < emi_mesh.l_numelem; eidx++)
  {
    // set element view to current element
    stiffness_view.set_elem(eidx);
 
    // calculate row/col indices of entries
    mesh_int_t nnodes = stiffness_view.num_nodes();
 
    row_idx.resize(nnodes*row_dpn);
    col_idx.resize(nnodes*col_dpn);
    canonic_indices<mesh_int_t,SF_int>(stiffness_view.nodes(), stiffness_petsc_nbr.data(), nnodes, row_dpn, row_idx.data());
    canonic_indices<mesh_int_t,SF_int>(stiffness_view.nodes(), stiffness_petsc_nbr.data(), nnodes, col_dpn, col_idx.data());
 
    // call integrator
    stiffness_integrator(stiffness_view, ebuff);
    ebuff *= stiffness_scale;
 
    // add values into system matrix
    const bool add = true;
    mat.set_values(row_idx, col_idx, ebuff.data(), add);
  }
 
  // Assemble the global surface mass matrix (M).
  // Each local mass matrix computed for a surface element is mapped into the global system
  // using the corresponding global row and column indices of that element’s nodes.
 
  const vector<SF_int> & ref_nbr = surface_mesh.get_numbering(SF::NBR_REF);
  const vector<SF_int> & petsc_nbr = surface_mesh.get_numbering(SF::NBR_PETSC);
 
  const SF::vector<mesh_int_t> & rnod_emi = emi_mesh.get_numbering(SF::NBR_REF);
  hashmap::unordered_map<T,T> g2l_emi;
  hashmap::unordered_map<T,T> l2g_emi;
  for(size_t i=0; i<rnod_emi.size(); i++){
    g2l_emi[rnod_emi[i]] = i;
    l2g_emi[i] = rnod_emi[i];
  }
  const SF::vector<mesh_int_t> & rnod = surface_mesh.get_numbering(SF::NBR_REF);
  hashmap::unordered_map<T,T> g2l;
  hashmap::unordered_map<T,T> l2g;
 
  for(size_t i=0; i<rnod.size(); i++){
    g2l[rnod[i]] = i;
    l2g[i] = rnod[i];
  }
 
  // start with assembly mass matrix into lhs matrix
  element_view<mesh_int_t, mesh_real_t> view(surface_mesh, NBR_PETSC);
  for(size_t eidx=0; eidx < surface_mesh.l_numelem; eidx++)
  {
    // set element view to current element
    view.set_elem(eidx);
    mesh_int_t tag = surface_mesh.tag[eidx]; // selected face
 
    // get the counter part
    std::vector<int> elem_nodes;  
    std::vector<int> elem_nodes_old;  
    std::vector<int> elem_nodes_first;
    std::vector<int> elem_nodes_second;
 
    std::vector<int> petsc_first;
    std::vector<int> petsc_second;
    
    T tag_first = 0;
    T tag_second = 0;
 
    for (int n = surface_mesh.dsp[eidx]; n < surface_mesh.dsp[eidx+1];n++)
    {
      T l_idx = surface_mesh.con[n];
 
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(l2g[l_idx],tag);
      mesh_int_t Index_new = map_vertex_tag_to_dof_petsc[Index_tag_old].first;
      elem_nodes.push_back(Index_new);
      elem_nodes_old.push_back(l2g[l_idx]);
    }
 
    std::sort(elem_nodes.begin(),elem_nodes.end());
    if(elem_nodes.size()==2){
      tuple_key key;
 
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      std::pair<tuple_value, tuple_value> value = line_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
    }
    else if(elem_nodes.size()==3){
      tri_key key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      std::pair<tri_value, tri_value> value = tri_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
    }
    else if(elem_nodes.size()==4){
      quad_key key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      key.v4 = elem_nodes[3];
      std::pair<quad_value, quad_value> value = quad_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
    }
 
    if(tag != tag_first)  std::swap(tag_first, tag_second);
    
    std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
 
    for (int indx = 0; indx < elem_nodes_old.size(); ++indx)
    {
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old_first;
      Index_tag_old_first = std::make_pair(elem_nodes_old[indx],tag_first);
      std::pair <mesh_int_t,mesh_int_t> newIndex_petsc_first = map_vertex_tag_to_dof_petsc[Index_tag_old_first];
      mesh_int_t newIndex_first = newIndex_petsc_first.first;
      petsc_first.push_back(newIndex_petsc_first.second);
 
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old_second;
      Index_tag_old_second = std::make_pair(elem_nodes_old[indx],tag_second);
      std::pair <mesh_int_t,mesh_int_t> newIndex_petsc_second = map_vertex_tag_to_dof_petsc[Index_tag_old_second];
      mesh_int_t newIndex_second = newIndex_petsc_second.first;
      petsc_second.push_back(newIndex_petsc_second.second);
    }
 
    // calculate row/col indices of entries
    SF_int nnodes = view.num_nodes();
    idx.resize(nnodes*row_dpn);
    idx_counter.resize(nnodes*col_dpn);
 
    {
      for(T i=0; i<nnodes; i++)
      {
        for(short j=0; j<row_dpn; j++){
          idx[i*row_dpn + j] = (petsc_first[i])*row_dpn + j;
        }
        for(short j=0; j<col_dpn; j++){
          idx_counter[i*col_dpn + j] = (petsc_second[i])*col_dpn + j;            
        }
      }
    }
 
    // add values into system matrix
    const bool add = true;
 
    // call integrator
    mass_integrator(view, ebuff);
    ebuff_mass.assign(nnodes, nnodes, 0.0);
    ebuff_mass = ebuff;
    // Note: surface mesh by default has both face and counter face; 
    // to avoid the dublication, factor 0.5 is necessary
    ebuff*=0.5 * mass_scale;
    ebuff_mass*=0.5;
                             
    mat.set_values(idx, idx, ebuff.data(), add);
    mat.set_values(idx_counter, idx_counter, ebuff.data(), add);
 
    mat_surf.set_values(idx, idx, ebuff_mass.data(), add);
    mat_surf.set_values(idx_counter, idx_counter, ebuff_mass.data(), add);
 
    ebuff*=(-1.0);
    ebuff_mass*=(-1.0);
 
    mat.set_values(idx, idx_counter, ebuff.data(), add);
    mat.set_values(idx_counter, idx, ebuff.data(), add);
 
    mat_surf.set_values(idx, idx_counter, ebuff_mass.data(), add);
    mat_surf.set_values(idx_counter, idx, ebuff_mass.data(), add);
 
    // Each face contributes the following 2x2 block structure:
    //
    //   s* [ + M    - M ;
    //        - M    + M ]
    //
    // with s = 0.5 * alpha,
    // where M is the local face mass matrix.
 
  }
  // finish assembly and progress output
  mat.finish_assembly();
  mat_surf.finish_assembly();
}
 
template<class tuple_key, class tuple_value, class tri_key, class tri_value, class quad_key, class quad_value, class T, class S>
inline void assign_resting_potential_from_ionic_models_on_myocyte(abstract_matrix<T,S> & Vm_emi,
                          abstract_vector<T, S>* vb,
                          SF::vector<mesh_int_t> elemTag_emi_mesh,
                          hashmap::unordered_map<std::pair<T,T>, 
                          std::pair<T,T>> map_vertex_tag_to_dof_petsc,
                          hashmap::unordered_map<tuple_key,
                          std::pair<tuple_value, tuple_value>> line_face,
                          hashmap::unordered_map<tri_key,
                          std::pair<tri_value, tri_value>> tri_face, 
                          hashmap::unordered_map<quad_key,
                          std::pair<quad_value, quad_value>> quad_face,
                          const meshdata<T,mesh_real_t> & surface_mesh,
                          const meshdata<T,mesh_real_t> & emi_mesh)
{
  // we want to make sure that the element integrator fits the matrix
  T row_dpn = 1; T col_dpn = 1;
 
  const vector<SF_int> & petsc_nbr = emi_mesh.get_numbering(NBR_PETSC);// remove it
  const SF::vector<T> & rnod_emi = emi_mesh.get_numbering(SF::NBR_REF);
  hashmap::unordered_map<T,T> g2l_emi;
  hashmap::unordered_map<T,T> l2g_emi;
  for(size_t i=0; i<rnod_emi.size(); i++){
    g2l_emi[rnod_emi[i]] = i;
    l2g_emi[i] = rnod_emi[i];
  }
 
  const SF::vector<T> & rnod = surface_mesh.get_numbering(SF::NBR_REF);
  const vector<SF_int> & petsc_nbr_surface = surface_mesh.get_numbering(NBR_PETSC);// remove it
 
  hashmap::unordered_map<T,T> g2l;
  hashmap::unordered_map<T,T> l2g;
 
  hashmap::unordered_map<T,SF_real> tag_2_vm;
 
  for(size_t i=0; i<rnod.size(); i++){
    g2l[rnod[i]] = i;
    l2g[i] = rnod[i];
  }
  auto vb_data = vb->const_ptr();
 
  // start with assembly
  for(size_t eidx=0; eidx < surface_mesh.l_numelem; eidx++)
  {
    T tag = surface_mesh.tag[eidx];
 
    // get the counter part
    std::vector<int> elem_nodes;
    bool tag_current;
    T tag_first = 0;
    T tag_second = 0;
    T mem_first = 0;
    T mem_second = 0;
 
    T sign = 1.0;
 
    for (int n = surface_mesh.dsp[eidx], i = 0; n < surface_mesh.dsp[eidx+1];n++,i++)
    {
      T l_idx = surface_mesh.con[n];
      std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
      Index_tag_old = std::make_pair(l2g[l_idx],tag);
      mesh_int_t Index_new = map_vertex_tag_to_dof_petsc[Index_tag_old].first;
      elem_nodes.push_back(Index_new);
    }
 
    std::sort(elem_nodes.begin(),elem_nodes.end());
    if(elem_nodes.size()==2){
      tuple_key key;
 
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      std::pair<tuple_value, tuple_value> value = line_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
 
      mem_first = value.first.mem;
      mem_second = value.second.mem;
    }
    else if(elem_nodes.size()==3){
      tri_key key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      std::pair<tri_value, tri_value> value = tri_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
 
      mem_first = value.first.mem;
      mem_second = value.second.mem;
    }
    else if(elem_nodes.size()==4){
      quad_key key;
      key.v1 = elem_nodes[0];
      key.v2 = elem_nodes[1];
      key.v3 = elem_nodes[2];
      key.v4 = elem_nodes[3];
      std::pair<quad_value, quad_value> value = quad_face[key];
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
 
      mem_first = value.first.mem;
      mem_second = value.second.mem;
    }
 
    SF_int nnodes = elem_nodes.size();
    for (int i = 0; i < nnodes; ++i)
    {
      if(mem_first==1 || mem_second==1){
        tag_2_vm[tag_first] = vb_data[eidx];
        tag_2_vm[tag_second] = vb_data[eidx];
      }
    }
  }
 
  vb->const_release_ptr(vb_data);
 
  // assign all the tags
  {
    // we want to make sure that the element integrator fits the matrix
    mesh_int_t row_dpn = 1, col_dpn = 1;
 
    // allocate row and value 
    vector<T> row_idx(SF_MAX_ELEM_NODES * row_dpn);
    vector<T> col_idx(SF_MAX_ELEM_NODES * row_dpn);
    vector<SF_real> val_idx(SF_MAX_ELEM_NODES * row_dpn);
 
    // allocate elem index buffer
    dmat<SF_real> ebuff(SF_MAX_ELEM_NODES * row_dpn, SF_MAX_ELEM_NODES * col_dpn);
 
    const vector<mesh_int_t> & petsc_nbr = emi_mesh.get_numbering(NBR_PETSC);
 
    // start with assembly
    element_view<mesh_int_t, mesh_real_t> view(emi_mesh, NBR_PETSC);
    for(size_t eidx=0; eidx < emi_mesh.l_numelem; eidx++)
    {
      // set element view to current element
      view.set_elem(eidx);
      T tag = emi_mesh.tag[eidx];
      // calculate row/col indices of entries
      mesh_int_t nnodes = view.num_nodes();
 
      row_idx.resize(nnodes*row_dpn);val_idx.resize(nnodes*row_dpn);
      col_idx.resize(nnodes*row_dpn);
      ebuff.assign(nnodes, nnodes, 0.0);
      canonic_indices<mesh_int_t,SF_int>(view.nodes(), petsc_nbr.data(), nnodes, row_dpn, row_idx.data());
      canonic_indices<mesh_int_t,SF_int>(view.nodes(), petsc_nbr.data(), nnodes, col_dpn, col_idx.data());
 
      for (int i = 0; i < nnodes; ++i)
      {
        ebuff[i][i] = 0;
        // elemTag_emi_mesh[eidx]==2 means that the tag belongs to the intracellular region
        if(elemTag_emi_mesh[eidx]==2){
          ebuff[i][i] = tag_2_vm[tag];
        }
      }
      Vm_emi.set_values(row_idx, col_idx, ebuff.data(), false);
    }
    Vm_emi.finish_assembly();
  }
}
 
}
 
#endif
#endif