sf_interface.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
// ----------------------------------------------------------------------------
 
#include "petsc_utils.h" // TODO: for EXIT
#include "sf_interface.h"
#include "sim_utils.h"
 
namespace opencarp {
 
sf_mesh & get_mesh(const mesh_t gt)
{
  auto it = user_globals::mesh_reg.find(gt);
  if(it == user_globals::mesh_reg.end()) {
    log_msg(0,5,0, "%s error: mesh of type \"%s\" not found! Aborting!",
            __func__, get_mesh_type_name(gt));
    EXIT(EXIT_FAILURE);
  }
 
  return it->second;
}
 
const char*
get_mesh_type_name(mesh_t t) {
  switch(t) {
    case intra_elec_msh: return "Intracellular Electric";
    case extra_elec_msh: return "Extracellular Electric";
    case eikonal_msh:    return "Eikonal";
    case elasticity_msh: return "Elasticity";
    case fluid_msh:      return "Fluid";
    case emi_msh:        return "EMI volume"; // includes both discontinuous faces and inner domain structures. 
    case emi_surface_msh: return "EMI surface"; // Interface mesh from EMI mesh (including gap junction face and membrane faces)
    case emi_surface_counter_msh: return "EMI surface with counter";  // The interface mesh is constructed from the EMI mesh, incorporating both gap junction faces and membrane faces. Additionally, each face includes its corresponding counter face between two subdomains, where current will be computed on both face and counter face 
    case emi_surface_unique_face_msh: return "EMI surface with single face on EMI interface"; // Interface mesh from EMI mesh (including gap junction face and membrane faces)
    case reference_msh:  return "Reference";
    case phie_recv_msh:  return "Phie Recovery";
    default: return NULL;
  }
}
 
bool mesh_is_registered(const mesh_t gt)
{
  return user_globals::mesh_reg.find(gt) != user_globals::mesh_reg.end();
}
 
SF::scattering*
register_scattering(const int from, const int to, const SF::SF_nbr nbr, const int dpn)
{
  int rank = get_rank();
 
  sf_mesh & from_mesh = get_mesh(mesh_t(from));
  // we compute the petsc numbering of only the algebraicly distributed nodes
  const SF::vector<mesh_int_t> & from_alg_nod   = from_mesh.pl.algebraic_nodes();
  const SF::vector<mesh_int_t> & from_nodal_nbr = from_mesh.get_numbering(nbr);
  SF::scattering* ret = NULL;
  SF::quadruple<int> spec = {from, to, nbr, dpn};
 
  if(to != ALG_TO_NODAL) {
    sf_mesh & to_mesh   = get_mesh(mesh_t(to));
    // we set up the inter-domain index mapping in the registry map_reg
    SF::index_mapping<mesh_int_t> & idx_map = user_globals::map_reg[spec];
    SF::inter_domain_mapping(from_mesh, to_mesh, nbr, idx_map);
 
    SF::vector<mesh_int_t> to_numbering(from_alg_nod.size());
    SF::vector<mesh_int_t> from_numbering(from_alg_nod.size());
    int err = 0;
 
    for(size_t i=0; i<from_alg_nod.size(); i++) {
      mesh_int_t from_idx = from_nodal_nbr[from_alg_nod[i]];
      mesh_int_t to_idx   = idx_map.forward_map(from_idx);
 
      if(to_idx < 0) {
        err++;
        break;
      }
 
      from_numbering[i] = from_idx;
      to_numbering[i] = to_idx;
    }
 
    if(get_global(err, MPI_SUM)) {
      log_msg(0,5,0, "%s error: Bad inter-domain mapping. Aborting!", __func__);
      EXIT(1);
    }
 
    ret = user_globals::scatter_reg.register_scattering(spec, from_mesh.pl.algebraic_layout(),
           to_mesh.pl.algebraic_layout(), from_numbering, to_numbering, rank, dpn);
 
  }
  else {
    const SF::vector<mesh_int_t> & nodal_layout = from_mesh.pl.layout();
    SF::vector<mesh_int_t>         to_numbering(from_nodal_nbr.size());
 
    SF::interval(to_numbering, nodal_layout[rank], nodal_layout[rank+1]);
 
    ret = user_globals::scatter_reg.register_scattering(spec, from_mesh.pl.algebraic_layout(),
           nodal_layout, from_nodal_nbr, to_numbering, rank, dpn);
  }
 
  return ret;
}
 
SF::scattering* get_scattering(const int from, const int to, const SF::SF_nbr nbr, const int dpn)
{
  SF::quadruple<int> spec = {from, to, int(nbr), dpn};
  SF::scattering* sc = user_globals::scatter_reg.get_scattering(spec);
 
  if(sc == NULL)
    return register_scattering(from, to, nbr, dpn);
  else
    return sc;
}
 
bool have_scattering(const int from, const int to, const SF::SF_nbr nbr, const int dpn)
{
  SF::quadruple<int> spec = {from, to, int(nbr), dpn};
  SF::scattering* sc = user_globals::scatter_reg.get_scattering(spec);
 
  return sc != NULL;
}
 
SF::scattering*
register_permutation(const int mesh_id, const int perm_id, const int dpn)
{
  sf_mesh & mesh = get_mesh(mesh_t(mesh_id));
  const SF::vector<mesh_int_t> & petsc_nbr = mesh.get_numbering(SF::NBR_PETSC);
  const SF::vector<mesh_int_t> & alg_nod   = mesh.pl.algebraic_nodes();
  SF::quadruple<int> spec = {mesh_id, perm_id, 0, dpn};
 
  switch(perm_id) {
    case PETSC_TO_CANONICAL:
    {
      const SF::vector<mesh_int_t> & canon_nbr = mesh.get_numbering(SF::NBR_SUBMESH);
      SF::vector<mesh_int_t> petsc_alg_nod(alg_nod.size()), canon_alg_nod(alg_nod.size());
 
      for(size_t i=0; i<alg_nod.size(); i++)
      {
        petsc_alg_nod[i] = petsc_nbr[alg_nod[i]];
        canon_alg_nod[i] = canon_nbr[alg_nod[i]];
      }
 
      return user_globals::scatter_reg.register_permutation(spec, petsc_alg_nod, canon_alg_nod,
             mesh.g_numpts, mesh.g_numpts, dpn);
    }
 
    case ELEM_PETSC_TO_CANONICAL:
    {
      const SF::vector<mesh_int_t> & canon_nbr = mesh.get_numbering(SF::NBR_ELEM_SUBMESH);
      const SF::vector<mesh_int_t> & elem_layout = mesh.epl.algebraic_layout();
      const int rank = get_rank();
 
      SF::vector<mesh_int_t> petsc_nod(canon_nbr.size()), canon_nod(canon_nbr.size());
 
      for(size_t i=0; i<canon_nbr.size(); i++)
      {
        petsc_nod[i] = elem_layout[rank] + i;
        canon_nod[i] = canon_nbr[i];
      }
 
      return user_globals::scatter_reg.register_permutation(spec, petsc_nod, canon_nod,
             mesh.g_numelem, mesh.g_numelem, dpn);
    }
/*
    case PETSC_TO_PT:
    {
      const SF::vector<mesh_int_t> & pt_nbr = mesh.get_numbering(SF::NBR_SOLVER);
      SF::vector<mesh_int_t> petsc_alg_nod(alg_nod.size()), pt_alg_nod(alg_nod.size());
 
      for(size_t i=0; i<alg_nod.size(); i++)
      {
        petsc_alg_nod[i] = petsc_nbr[alg_nod[i]];
        pt_alg_nod[i]    = pt_nbr[alg_nod[i]];
      }
 
      return user_globals::scatter_reg.register_permutation(spec, petsc_alg_nod, pt_alg_nod,
             mesh.g_numpts, mesh.g_numpts, dpn);
    }
*/
    default:
      log_msg(0,5,0, "%s error: Unknown permutation id. Aborting!", __func__);
      EXIT(1);
  }
}
 
SF::scattering*
get_permutation(const int mesh_id, const int perm_id, const int dpn)
{
  SF::quadruple<int> spec = {mesh_id, perm_id, 0, dpn};
  SF::scattering* sc = user_globals::scatter_reg.get_scattering(spec);
 
  if(sc == NULL)
    return register_permutation(mesh_id, perm_id, dpn);
  else
    return sc;
}
 
bool have_permutation(const int mesh_id, const int perm_id, const int dpn)
{
  SF::quadruple<int> spec = {mesh_id, perm_id, 0, dpn};
  SF::scattering* sc = user_globals::scatter_reg.get_scattering(spec);
 
  return sc != NULL;
}
 
int get_phys_index(int physreg)
{
  int idx = -1;
  for(int i=0; i<param_globals::num_phys_regions; i++) {
    if(param_globals::phys_region[i].ptype == physreg) {
      idx = i;
      break;
    }
  }
  return idx;
}
 
bool phys_defined(int physreg)
{
  int idx = get_phys_index(physreg);
  return idx > -1;
}
 
 
}  // namespace opencarp