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 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;
}

}  // namespace opencarp