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"
#include "SF_init.h"
#include "mpi_utils.h"
#include <array>
#include <cmath>
#include <limits>
#include <stdexcept>
#include <type_traits>
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, std::array<double, SF_MAX_ELEM_NODES>& interpolationCoefficients, double & length)
{
  SF::Point d = p2 - p1;
  double lengthSquared = inner_prod(d, d);
  length = std::sqrt(lengthSquared);
 
  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, std::array<double, SF_MAX_ELEM_NODES>& 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, std::array<double, SF_MAX_ELEM_NODES>& 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;
 
  std::array<double, SF_MAX_ELEM_NODES> interpolationCoefficients{};
  double 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);
  }
 
  // Use the same face shape weights for the scatter mass as for interpolation.
  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*area*interpolationCoefficients[0];
    ebuff_s[0][1] = mass_scale*area*interpolationCoefficients[1];
  }
  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*area*interpolationCoefficients[0];
    ebuff_s[0][1] = mass_scale*area*interpolationCoefficients[1];
    ebuff_s[0][2] = mass_scale*area*interpolationCoefficients[2];
  }
  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*area*interpolationCoefficients[0];
    ebuff_s[0][1] = mass_scale*area*interpolationCoefficients[1];
    ebuff_s[0][2] = mass_scale*area*interpolationCoefficients[2];
    ebuff_s[0][3] = mass_scale*area*interpolationCoefficients[3];
  }
}
 
 
template<class T, class S, class emi_index_rank>
inline void construct_direct_unique_both_operators(
    SF::abstract_matrix<T, S>*& operator_unique_to_both_faces,
    SF::abstract_matrix<T, S>*& operator_both_to_unique_face,
    hashmap::unordered_map<mesh_int_t, std::pair<emi_index_rank, emi_index_rank>>& map_elem_uniqueFace_to_elem_bothface,
    const hashmap::unordered_map<mesh_int_t, std::pair<emi_index_rank, emi_index_rank>>& map_elem_uniqueFace_to_elem_oneface,
    const SF::vector<mesh_int_t>& vec_both_to_one_face,
    const SF::meshdata<mesh_int_t, mesh_real_t>& emi_surfmesh_w_counter_face,
    const SF::meshdata<mesh_int_t, mesh_real_t>& emi_surfmesh_unique_face,
    int max_row_entries_emi,
    int dpn,
    typename SF::abstract_vector<T, S>::ltype alg_surface_type)
{
  // Global/local dimensions for the direct unique-face <-> both-face operators.
  T M = emi_surfmesh_w_counter_face.g_numelem;
  T m = emi_surfmesh_w_counter_face.l_numelem;
  T M_unique_face = emi_surfmesh_unique_face.g_numelem;
  T m_unique_face = emi_surfmesh_unique_face.l_numelem;
 
  int rank = 0;
  int comm_size = 0;
  MPI_Comm_rank(emi_surfmesh_w_counter_face.comm, &rank);
  MPI_Comm_size(emi_surfmesh_w_counter_face.comm, &comm_size);
 
  SF::vector<long int> layout;
  SF::layout_from_count<long int>(emi_surfmesh_w_counter_face.l_numelem, layout, emi_surfmesh_w_counter_face.comm);
  T m_l = layout[rank];
 
  SF::vector<long int> layout_unique_face;
  SF::layout_from_count<long int>(emi_surfmesh_unique_face.l_numelem, layout_unique_face, emi_surfmesh_unique_face.comm);
  T m_unique_face_l = layout_unique_face[rank];
 
  // Rebuild the direct mapping from scratch every time the operators are assembled.
  map_elem_uniqueFace_to_elem_bothface.clear();
 
  // Gather the local both->one lookup so every rank can recover the matching
  // both-face indices for remote one-face owners.
  std::vector<int> both_counts(comm_size, 0), both_displs(comm_size, 0);
  int local_both_count = static_cast<int>(vec_both_to_one_face.size());
  MPI_Allgather(&local_both_count, 1, MPI_INT, both_counts.data(), 1, MPI_INT,
                emi_surfmesh_w_counter_face.comm);
 
  int total_both_count = 0;
  for (int i = 0; i < comm_size; i++) {
    both_displs[i] = total_both_count;
    total_both_count += both_counts[i];
  }
 
  std::vector<mesh_int_t> all_both_to_one(total_both_count);
  // Exchange raw bytes so the gather stays correct even if mesh_int_t differs
  // from int on a given build.
  std::vector<int> both_byte_counts(comm_size, 0), both_byte_displs(comm_size, 0);
  for (int i = 0; i < comm_size; i++) {
    both_byte_counts[i] = both_counts[i] * static_cast<int>(sizeof(mesh_int_t));
    both_byte_displs[i] = both_displs[i] * static_cast<int>(sizeof(mesh_int_t));
  }
  const int local_both_bytes = local_both_count * static_cast<int>(sizeof(mesh_int_t));
  MPI_Allgatherv(reinterpret_cast<const unsigned char*>(vec_both_to_one_face.data()),
                 local_both_bytes, MPI_BYTE,
                 reinterpret_cast<unsigned char*>(all_both_to_one.data()),
                 both_byte_counts.data(), both_byte_displs.data(), MPI_BYTE,
                 emi_surfmesh_w_counter_face.comm);
 
  // Build per-rank one->both lookup tables. A one-face entry can map to up to
  // two both-face entries, so we store first and second occurrences explicitly.
  std::vector<std::vector<mesh_int_t>> one_to_both_first(comm_size);
  std::vector<std::vector<mesh_int_t>> one_to_both_second(comm_size);
  for (int r = 0; r < comm_size; r++) {
    mesh_int_t max_one = -1;
    for (int i = 0; i < both_counts[r]; i++) {
      mesh_int_t one_idx = all_both_to_one[both_displs[r] + i];
      if (one_idx > max_one) max_one = one_idx;
    }
    if (max_one < 0) continue;
 
    one_to_both_first[r].assign(max_one + 1, -1);
    one_to_both_second[r].assign(max_one + 1, -1);
    for (int i = 0; i < both_counts[r]; i++) {
      mesh_int_t one_idx = all_both_to_one[both_displs[r] + i];
      if (one_idx < 0) continue;
      if (one_to_both_first[r][one_idx] < 0) one_to_both_first[r][one_idx] = i;
      else one_to_both_second[r][one_idx] = i;
    }
  }
 
  // Convert the stored unique->one relation into a direct unique->both relation.
  for (const auto& [unique_idx, one_face_pair] : map_elem_uniqueFace_to_elem_oneface) {
    const auto& first_oneface = one_face_pair.first;
    const auto& second_oneface = one_face_pair.second;
 
    auto map_one_to_both = [&](const emi_index_rank& one_face_idx, bool use_second) {
      emi_index_rank both_face_idx;
      both_face_idx.index = -1;
      both_face_idx.rank = -1;
 
      if (one_face_idx.rank < 0 || one_face_idx.rank >= comm_size || one_face_idx.index < 0) {
        return both_face_idx;
      }
      if (one_face_idx.index >= static_cast<int>(one_to_both_first[one_face_idx.rank].size())) {
        return both_face_idx;
      }
 
      mesh_int_t both_idx = one_to_both_first[one_face_idx.rank][one_face_idx.index];
      if (use_second &&
          one_face_idx.index < static_cast<int>(one_to_both_second[one_face_idx.rank].size()) &&
          one_to_both_second[one_face_idx.rank][one_face_idx.index] >= 0) {
        both_idx = one_to_both_second[one_face_idx.rank][one_face_idx.index];
      }
 
      if (both_idx >= 0) {
        both_face_idx.index = both_idx;
        both_face_idx.rank = one_face_idx.rank;
      }
      return both_face_idx;
    };
 
    const bool same_oneface =
        (first_oneface.index >= 0 && second_oneface.index >= 0 &&
         first_oneface.index == second_oneface.index &&
         first_oneface.rank == second_oneface.rank);
 
    map_elem_uniqueFace_to_elem_bothface[unique_idx] = std::make_pair(
        map_one_to_both(first_oneface, false),
        map_one_to_both(second_oneface, same_oneface));
  }
 
  // Assemble the expansion operator unique -> both. Each unique-face column
  // writes into one or two both-face rows, depending on ownership/splitting.
  SF::init_matrix(&operator_unique_to_both_faces);
  operator_unique_to_both_faces->init(M, M_unique_face, m, m_unique_face, m_l, 1);
  operator_unique_to_both_faces->zero();
 
  // Broadcast owner-provided row locations once so each rank can insert its
  // owned rows without repeating the communication during exact-preallocation
  // replay. This is a known remaining EMI scalability limit for 64-bit builds:
  // see docs/64BIT_SUPPORT.md.
  if (emi_surfmesh_unique_face.g_numelem > static_cast<size_t>(std::numeric_limits<int>::max())) {
    throw std::runtime_error(
        "EMI unique-face transfer operator currently requires fewer than INT_MAX unique faces");
  }
  const int M_unique_global = static_cast<int>(emi_surfmesh_unique_face.g_numelem);
  std::vector<int> first_rank(M_unique_global, -1), first_idx(M_unique_global, -1);
  std::vector<int> second_rank(M_unique_global, -1), second_idx(M_unique_global, -1);
 
  for (const auto& [local_unique_idx, both_pair] : map_elem_uniqueFace_to_elem_bothface) {
    int global_unique = static_cast<int>(layout_unique_face[rank] + local_unique_idx);
    if (global_unique < 0 || global_unique >= M_unique_global) continue;
 
    first_rank[global_unique] = static_cast<int>(both_pair.first.rank);
    first_idx[global_unique] = static_cast<int>(both_pair.first.index);
    second_rank[global_unique] = static_cast<int>(both_pair.second.rank);
    second_idx[global_unique] = static_cast<int>(both_pair.second.index);
  }
 
  MPI_Allreduce(MPI_IN_PLACE, first_rank.data(), M_unique_global, MPI_INT, MPI_MAX, emi_surfmesh_w_counter_face.comm);
  MPI_Allreduce(MPI_IN_PLACE, first_idx.data(), M_unique_global, MPI_INT, MPI_MAX, emi_surfmesh_w_counter_face.comm);
  MPI_Allreduce(MPI_IN_PLACE, second_rank.data(), M_unique_global, MPI_INT, MPI_MAX, emi_surfmesh_w_counter_face.comm);
  MPI_Allreduce(MPI_IN_PLACE, second_idx.data(), M_unique_global, MPI_INT, MPI_MAX, emi_surfmesh_w_counter_face.comm);
 
  auto assemble_unique_to_both = [&]() {
    SF::vector<SF_int> row_idx(1), col_idx(1);
    SF::dmat<SF_real> ebuff(1, 1);
    ebuff.assign(1, 1, 1.0);
 
    for (int gid = 0; gid < M_unique_global; gid++) {
      col_idx[0] = gid;
      if (first_rank[gid] == rank && first_idx[gid] >= 0) {
        row_idx[0] = static_cast<SF_int>(layout[rank] + first_idx[gid]);
        operator_unique_to_both_faces->set_values(row_idx, col_idx, ebuff.data(), false);
      }
      if (second_rank[gid] == rank && second_idx[gid] >= 0) {
        row_idx[0] = static_cast<SF_int>(layout[rank] + second_idx[gid]);
        operator_unique_to_both_faces->set_values(row_idx, col_idx, ebuff.data(), false);
      }
    }
    operator_unique_to_both_faces->finish_assembly();
  };
 
  if (operator_unique_to_both_faces->begin_exact_preallocation()) {
    assemble_unique_to_both();
    operator_unique_to_both_faces->finalize_exact_preallocation();
  }
  assemble_unique_to_both();
 
  // Assemble the restriction/averaging operator both -> unique from the direct map.
  SF::init_matrix(&operator_both_to_unique_face);
  operator_both_to_unique_face->init(M_unique_face, M, m_unique_face, m, m_unique_face_l, 2);
  operator_both_to_unique_face->zero();
  if (operator_both_to_unique_face->begin_exact_preallocation()) {
    assemble_map_both_to_unique(*operator_both_to_unique_face,
                                map_elem_uniqueFace_to_elem_bothface,
                                emi_surfmesh_unique_face,
                                emi_surfmesh_w_counter_face);
    operator_both_to_unique_face->finalize_exact_preallocation();
  }
  assemble_map_both_to_unique(*operator_both_to_unique_face,
                                  map_elem_uniqueFace_to_elem_bothface,
                                  emi_surfmesh_unique_face,
                                  emi_surfmesh_w_counter_face);
 
#ifdef EMI_DEBUG_MESH
  {
    // Validate mapping coverage and check whether the direct operators compose
    // to the identity on a constant unique-face field.
    mesh_int_t local_valid_first = 0, local_valid_second = 0;
    mesh_int_t local_same_column = 0;
    SF::vector<long int> layout_both_dbg;
    SF::layout_from_count<long int>(emi_surfmesh_w_counter_face.l_numelem, layout_both_dbg, emi_surfmesh_w_counter_face.comm);
    SF::vector<long int> layout_unique_dbg;
    SF::layout_from_count<long int>(emi_surfmesh_unique_face.l_numelem, layout_unique_dbg, emi_surfmesh_unique_face.comm);
 
    for (const auto& [uidx, both_pair] : map_elem_uniqueFace_to_elem_bothface) {
      if (both_pair.first.index >= 0) local_valid_first++;
      if (both_pair.second.index >= 0) local_valid_second++;
 
      if (both_pair.first.index >= 0 && both_pair.second.index >= 0 &&
          both_pair.first.rank >= 0 && both_pair.second.rank >= 0) {
        mesh_int_t global_first = layout_both_dbg[both_pair.first.rank] + both_pair.first.index;
        mesh_int_t global_second = layout_both_dbg[both_pair.second.rank] + both_pair.second.index;
        if (global_first == global_second) local_same_column++;
      }
    }
    mesh_int_t global_valid_first = local_valid_first;
    mesh_int_t global_valid_second = local_valid_second;
    mesh_int_t global_same_column = local_same_column;
    const MPI_Datatype mesh_mpi_t = opencarp::mpi_datatype<mesh_int_t>();
    MPI_Allreduce(MPI_IN_PLACE, &global_valid_first, 1, mesh_mpi_t, MPI_SUM, emi_surfmesh_w_counter_face.comm);
    MPI_Allreduce(MPI_IN_PLACE, &global_valid_second, 1, mesh_mpi_t, MPI_SUM, emi_surfmesh_w_counter_face.comm);
    MPI_Allreduce(MPI_IN_PLACE, &global_same_column, 1, mesh_mpi_t, MPI_SUM, emi_surfmesh_w_counter_face.comm);
 
    int dbg_rank = -1;
    MPI_Comm_rank(emi_surfmesh_w_counter_face.comm, &dbg_rank);
    if (dbg_rank == 0) {
      log_msg(NULL, 0, 0, "DEBUG map_unique_to_both: valid_first=%jd valid_second=%jd (global M=%zu, M_unique=%zu)",
              opencarp::printable_int(global_valid_first),
              opencarp::printable_int(global_valid_second),
              (size_t)emi_surfmesh_w_counter_face.g_numelem,
              (size_t)emi_surfmesh_unique_face.g_numelem);
      log_msg(NULL, 0, 0, "DEBUG map_unique_to_both: same_column=%jd", opencarp::printable_int(global_same_column));
    }
 
    {
      std::vector<char> present(emi_surfmesh_unique_face.l_numelem, 0);
      for (const auto& [uidx, both_pair] : map_elem_uniqueFace_to_elem_bothface) {
        if (uidx >= 0 && uidx < (mesh_int_t)present.size()) present[uidx] = 1;
      }
      std::vector<int> missing;
      for (mesh_int_t i = 0; i < (mesh_int_t)present.size(); i++) {
        if (!present[i]) missing.push_back((int)i);
      }
      int local_missing = (int)missing.size();
      int comm_size_dbg = 0;
      MPI_Comm_size(emi_surfmesh_w_counter_face.comm, &comm_size_dbg);
      std::vector<int> counts(comm_size_dbg, 0), displs(comm_size_dbg, 0);
      MPI_Gather(&local_missing, 1, MPI_INT,
                 dbg_rank == 0 ? counts.data() : nullptr, 1, MPI_INT,
                 0, emi_surfmesh_w_counter_face.comm);
      if (dbg_rank == 0) {
        int total = 0;
        for (int r = 0; r < comm_size_dbg; r++) {
          displs[r] = total;
          total += counts[r];
        }
        std::vector<int> all_missing(total, -1);
        MPI_Gatherv(missing.data(), local_missing, MPI_INT,
                    all_missing.data(), counts.data(), displs.data(), MPI_INT,
                    0, emi_surfmesh_w_counter_face.comm);
        int offset = 0;
        bool any = false;
        for (int r = 0; r < comm_size_dbg; r++) {
          if (counts[r] == 0) {
            offset += counts[r];
            continue;
          }
          any = true;
          log_msg(NULL, 0, 0, "DEBUG unique missing: rank=%d count=%d", r, counts[r]);
          int to_print = counts[r] < 10 ? counts[r] : 10;
          for (int i = 0; i < to_print; i++) {
            log_msg(NULL, 0, 0, "DEBUG unique missing: rank=%d local_unique=%d", r, all_missing[offset + i]);
          }
          offset += counts[r];
        }
        if (!any) log_msg(NULL, 0, 0, "DEBUG unique missing: none");
      } else {
        MPI_Gatherv(missing.data(), local_missing, MPI_INT,
                    nullptr, nullptr, nullptr, MPI_INT,
                    0, emi_surfmesh_w_counter_face.comm);
      }
    }
 
    SF::abstract_vector<T, S>* dbg_unique = nullptr;
    SF::abstract_vector<T, S>* dbg_both = nullptr;
    SF::abstract_vector<T, S>* dbg_back = nullptr;
    SF::init_vector(&dbg_unique, emi_surfmesh_unique_face, dpn, alg_surface_type);
    SF::init_vector(&dbg_both, emi_surfmesh_w_counter_face, dpn, alg_surface_type);
    SF::init_vector(&dbg_back, emi_surfmesh_unique_face, dpn, alg_surface_type);
    dbg_unique->set(1.0);
    operator_unique_to_both_faces->mult(*dbg_unique, *dbg_both);
    operator_both_to_unique_face->mult(*dbg_both, *dbg_back);
 
    SF_real local_max_err = 0.0;
    mesh_int_t local_half_count = 0;
    SF_real* p = dbg_back->ptr();
    for (mesh_int_t i = 0; i < dbg_back->lsize(); i++) {
      SF_real err = std::abs(p[i] - 1.0);
      if (err > local_max_err) local_max_err = err;
      if (std::abs(p[i] - 0.5) < 1e-12) local_half_count++;
    }
    dbg_back->release_ptr(p);
    SF_real global_max_err = 0.0;
    mesh_int_t global_half_count = 0;
    MPI_Allreduce(&local_max_err, &global_max_err, 1, opencarp::mpi_datatype<SF_real>(),
                  MPI_MAX, emi_surfmesh_w_counter_face.comm);
    global_half_count = local_half_count;
    MPI_Allreduce(MPI_IN_PLACE, &global_half_count, 1, mesh_mpi_t, MPI_SUM,
                  emi_surfmesh_w_counter_face.comm);
    if (dbg_rank == 0) {
      log_msg(NULL, 0, 0, "DEBUG map_unique_to_both consistency: max_err=%.6e half_count=%jd",
              (double)global_max_err, opencarp::printable_int(global_half_count));
    }
 
    {
      const int max_print = 5;
      const int fields = 9;
      std::vector<int> local_buf(max_print * fields, -1);
      int filled = 0;
 
      for (const auto& [local_unique_idx, both_pair] : map_elem_uniqueFace_to_elem_bothface) {
        if (filled >= max_print) break;
        const auto& first_both = both_pair.first;
        const auto& second_both = both_pair.second;
 
        const bool first_valid = (first_both.index >= 0 && first_both.rank >= 0);
        const bool second_valid = (second_both.index >= 0 && second_both.rank >= 0);
        int count = (first_valid ? 1 : 0) + (second_valid ? 1 : 0);
        if (count != 1) continue;
 
        mesh_int_t global_unique = layout_unique_dbg[dbg_rank] + local_unique_idx;
        mesh_int_t global_first = first_valid ? (layout_both_dbg[first_both.rank] + first_both.index) : -1;
        mesh_int_t global_second = second_valid ? (layout_both_dbg[second_both.rank] + second_both.index) : -1;
 
        int base = filled * fields;
        local_buf[base + 0] = dbg_rank;
        local_buf[base + 1] = (int)local_unique_idx;
        local_buf[base + 2] = (int)global_unique;
        local_buf[base + 3] = (int)first_both.rank;
        local_buf[base + 4] = (int)first_both.index;
        local_buf[base + 5] = (int)global_first;
        local_buf[base + 6] = (int)second_both.rank;
        local_buf[base + 7] = (int)second_both.index;
        local_buf[base + 8] = (int)global_second;
        filled++;
      }
 
      int comm_size = 0;
      MPI_Comm_size(emi_surfmesh_w_counter_face.comm, &comm_size);
      std::vector<int> all_buf;
      if (dbg_rank == 0) all_buf.resize(comm_size * max_print * fields, -1);
 
      MPI_Gather(local_buf.data(), max_print * fields, MPI_INT,
                 dbg_rank == 0 ? all_buf.data() : nullptr, max_print * fields, MPI_INT,
                 0, emi_surfmesh_w_counter_face.comm);
 
      if (dbg_rank == 0) {
        for (int r = 0; r < comm_size; r++) {
          for (int i = 0; i < max_print; i++) {
            int base = (r * max_print + i) * fields;
            if (all_buf[base + 1] < 0) continue;
            log_msg(NULL, 0, 0,
                    "DEBUG map_unique_to_both bad: rank=%d local_unique=%d global_unique=%d "
                    "first(rank=%d idx=%d glob=%d) second(rank=%d idx=%d glob=%d)",
                    all_buf[base + 0], all_buf[base + 1], all_buf[base + 2],
                    all_buf[base + 3], all_buf[base + 4], all_buf[base + 5],
                    all_buf[base + 6], all_buf[base + 7], all_buf[base + 8]);
          }
        }
      }
    }
 
    delete dbg_unique;
    delete dbg_both;
    delete dbg_back;
  }
#endif
}
 
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,
                                        const SF::vector<mesh_int_t> & elemTag_surface_mesh,
                                        const 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,
                                        const hashmap::unordered_map<tuple_key,
                                        std::pair<tuple_value, tuple_value>> & line_face,
                                        const hashmap::unordered_map<tri_key,
                                        std::pair<tri_value, tri_value>> & tri_face,
                                        const 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 SF::vector<mesh_int_t> & rnod_emi = emi_mesh.get_numbering(SF::NBR_REF);
  hashmap::unordered_map<T,T> g2l_emi;
  g2l_emi.reserve(rnod_emi.size());
  for(size_t i=0; i<rnod_emi.size(); i++){
    g2l_emi[rnod_emi[i]] = i;
  }
 
 
  const SF::vector<mesh_int_t> & rnod = surface_mesh.get_numbering(SF::NBR_REF);
 
  auto find_petsc_index = [&](const std::pair<mesh_int_t,mesh_int_t>& key) -> const std::pair<mesh_int_t,mesh_int_t>& {
    auto it = map_vertex_tag_to_dof_petsc.find(key);
    if (it == map_vertex_tag_to_dof_petsc.end()) {
      std::cerr << "ERROR: PETSc index not found for vertex=" << key.first
                << " tag=" << key.second << std::endl;
      throw std::runtime_error("PETSc index not found for vertex/tag pair");
    }
    return it->second;
  };
 
  // 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);
 
  std::vector<mesh_int_t> elem_nodes;
  std::vector<mesh_int_t> elem_nodes_old;
  std::vector<mesh_int_t> petsc_first;
  std::vector<mesh_int_t> petsc_second;
  vector<SF::Point> face_coordinates;
  elem_nodes.reserve(SF_MAX_ELEM_NODES);
  elem_nodes_old.reserve(SF_MAX_ELEM_NODES);
  petsc_first.reserve(SF_MAX_ELEM_NODES);
  petsc_second.reserve(SF_MAX_ELEM_NODES);
  face_coordinates.reserve(SF_MAX_ELEM_NODES);
 
  // start with assembly
  for(size_t eidx=0; eidx < surface_mesh.l_numelem; eidx++)
  {
    T tag = surface_mesh.tag[eidx];
 
    // get the counter part
    elem_nodes.clear();
    elem_nodes_old.clear();
    petsc_first.clear();
    petsc_second.clear();
 
    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]; 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(rnod[l_idx],tag);
      mesh_int_t Index_new = find_petsc_index(Index_tag_old).first;
      elem_nodes.push_back(Index_new);
      elem_nodes_old.push_back(rnod[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];
      auto it = line_face.find(key);
      if (it == line_face.end()) throw std::runtime_error("Line interface face not found");
      const std::pair<tuple_value, tuple_value> & value = it->second;
 
      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];
      auto it = tri_face.find(key);
      if (it == tri_face.end()) throw std::runtime_error("Triangular interface face not found");
      const std::pair<tri_value, tri_value> & value = it->second;
 
      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];
      auto it = quad_face.find(key);
      if (it == quad_face.end()) throw std::runtime_error("Quadrilateral interface face not found");
      const std::pair<quad_value, quad_value> & value = it->second;
 
      tag_first = value.first.tag;
      tag_second = value.second.tag;
 
      mem_first = value.first.mem;
      mem_second = value.second.mem;
 
    }
 
    if(tag != tag_first) {
      std::swap(tag_first, tag_second);
      std::swap(mem_first, mem_second);
    }
 
    std::pair <mesh_int_t,mesh_int_t> Index_tag_old;
 
    for (size_t 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();
    face_coordinates.resize(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(rnod[l_idx],tag);
      mesh_int_t Index_new = find_petsc_index(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
 
    // BsM is a volume-row / face-column scatter. Insert entries explicitly
    // because row_idx contains the same face row repeated for B/Bi assembly.
    ebuff_s*=(sign);
    for(SF_int i = 0; i < nnodes; i++) {
      BsM.set_value(col_idx[i], row_idx[0], ebuff_s[0][i], 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,
                             const 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,
                             const hashmap::unordered_map<tuple_key,
                             std::pair<tuple_value, tuple_value>> & line_face,
                             const hashmap::unordered_map<tri_key,
                             std::pair<tri_value, tri_value>> & tri_face,
                             const 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 SF::vector<mesh_int_t> & rnod = surface_mesh.get_numbering(SF::NBR_REF);
 
  auto find_petsc_index = [&](const std::pair<mesh_int_t,mesh_int_t>& key) -> const std::pair<mesh_int_t,mesh_int_t>& {
    auto it = map_vertex_tag_to_dof_petsc.find(key);
    if (it == map_vertex_tag_to_dof_petsc.end()) {
      std::cerr << "ERROR: PETSc index not found for vertex=" << key.first
                << " tag=" << key.second << std::endl;
      throw std::runtime_error("PETSc index not found for vertex/tag pair");
    }
    return it->second;
  };
 
  // 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<mesh_int_t> elem_nodes;
    std::vector<mesh_int_t> elem_nodes_old;
    std::vector<mesh_int_t> elem_nodes_first;
    std::vector<mesh_int_t> elem_nodes_second;
 
    std::vector<mesh_int_t> petsc_first;
    std::vector<mesh_int_t> 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(rnod[l_idx],tag);
      mesh_int_t Index_new = find_petsc_index(Index_tag_old).first;
      elem_nodes.push_back(Index_new);
      elem_nodes_old.push_back(rnod[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];
      auto it = line_face.find(key);
      if (it == line_face.end()) throw std::runtime_error("Line interface face not found");
      const std::pair<tuple_value, tuple_value> & value = it->second;
 
      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];
      auto it = tri_face.find(key);
      if (it == tri_face.end()) throw std::runtime_error("Triangular interface face not found");
      const std::pair<tri_value, tri_value> & value = it->second;
 
      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];
      auto it = quad_face.find(key);
      if (it == quad_face.end()) throw std::runtime_error("Quadrilateral interface face not found");
      const std::pair<quad_value, quad_value> & value = it->second;
 
      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 (size_t 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 = find_petsc_index(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 = find_petsc_index(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_vector<T,S> & ui,
                          abstract_vector<T, S>* vb,
                          const SF::vector<mesh_int_t> & elemTag_emi_mesh,
                          const hashmap::unordered_map<std::pair<T,T>,
                          std::pair<T,T>> & map_vertex_tag_to_dof_petsc,
                          const hashmap::unordered_map<tuple_key,
                          std::pair<tuple_value, tuple_value>> & line_face,
                          const hashmap::unordered_map<tri_key,
                          std::pair<tri_value, tri_value>> & tri_face,
                          const 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)
{
  T row_dpn = 1;
 
  const SF::vector<T> & rnod = surface_mesh.get_numbering(SF::NBR_REF);
 
  hashmap::unordered_map<T,SF_real> tag_2_vm;
 
  auto find_petsc_index = [&](const std::pair<T,T>& key) -> const std::pair<T,T>& {
    auto it = map_vertex_tag_to_dof_petsc.find(key);
    if (it == map_vertex_tag_to_dof_petsc.end()) {
      std::cerr << "ERROR: PETSc index not found for vertex=" << key.first
                << " tag=" << key.second << std::endl;
      throw std::runtime_error("PETSc index not found for vertex/tag pair");
    }
    return it->second;
  };
 
  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<mesh_int_t> elem_nodes;
    T tag_first = 0;
    T tag_second = 0;
    T mem_first = 0;
    T mem_second = 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 <T,T> Index_tag_old;
      Index_tag_old = std::make_pair(rnod[l_idx],tag);
      mesh_int_t Index_new = find_petsc_index(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];
      auto it = line_face.find(key);
      if (it == line_face.end()) throw std::runtime_error("Line interface face not found");
      const std::pair<tuple_value, tuple_value> & value = it->second;
 
      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];
      auto it = tri_face.find(key);
      if (it == tri_face.end()) throw std::runtime_error("Triangular interface face not found");
      const std::pair<tri_value, tri_value> & value = it->second;
 
      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];
      auto it = quad_face.find(key);
      if (it == quad_face.end()) throw std::runtime_error("Quadrilateral interface face not found");
      const std::pair<quad_value, quad_value> & value = it->second;
 
      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 resting potentials directly into the owned ui entries. This replaces
  // the old diagonal Vm_myocyte_emi matrix, whose diagonal was immediately
  // copied into ui and never used as a real operator.
  {
    vector<T> row_idx(SF_MAX_ELEM_NODES * row_dpn);
    const vector<mesh_int_t> & petsc_nbr = emi_mesh.get_numbering(NBR_PETSC);
    T start = 0, stop = 0;
    ui.get_ownership_range(start, stop);
    ui.set(0.0);
    S* ui_data = ui.ptr();
 
    element_view<mesh_int_t, mesh_real_t> view(emi_mesh, NBR_PETSC);
    for(size_t eidx=0; eidx < emi_mesh.l_numelem; eidx++)
    {
      view.set_elem(eidx);
      T tag = emi_mesh.tag[eidx];
      mesh_int_t nnodes = view.num_nodes();
 
      row_idx.resize(nnodes*row_dpn);
      canonic_indices<mesh_int_t,SF_int>(view.nodes(), petsc_nbr.data(), nnodes, row_dpn, row_idx.data());
 
      for (int i = 0; i < nnodes; ++i)
      {
        // elemTag_emi_mesh[eidx]==2 means that the tag belongs to the intracellular region
        if(elemTag_emi_mesh[eidx]==2){
          if(row_idx[i] >= start && row_idx[i] < stop) {
            auto vm_it = tag_2_vm.find(tag);
            if (vm_it != tag_2_vm.end()) ui_data[row_idx[i] - start] = vm_it->second;
          }
        }
      }
    }
    ui.release_ptr(ui_data);
    ui.finish_assembly();
  }
}
 
}
 
#endif
#endif