SF_mesh_io_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_MESH_IO_EMI
#define _SF_MESH_IO_EMI
 
#include "SF_mesh_io.h"
#include "petsc_utils.h"
 
namespace SF {
 
template<class T, class S>
inline void insert_points_ptsData_to_dof( const meshdata<T, S> & mesh_original, meshdata<T, S> & mesh, 
                                          const SF_nbr numbering, 
                                          hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex,
                                          hashmap::unordered_map<T,T> & vertex2ptsdata,
                                          hashmap::unordered_map<T,T> & dof2ptsData)
{
  MPI_Comm comm = mesh_original.comm;
 
  const SF::vector<mesh_int_t> & rnod_original = mesh_original.get_numbering(numbering);
  hashmap::unordered_map<T,T> g2l_orig;
  for(size_t i=0; i<rnod_original.size(); i++){
    g2l_orig[rnod_original[i]] = i;
  }
 
  const SF::vector<mesh_int_t> & rnod = mesh.get_numbering(numbering);
  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];
  }
 
  mesh_int_t gmax = global_max(rnod, comm);
  mesh.g_numpts = gmax+1;
 
  SF::vector<mesh_real_t> Buffxyz_orig = mesh_original.xyz;
 
  mesh.xyz.resize(mesh.l_numpts*3);
 
  SF::meshdata<mesh_int_t, mesh_real_t>     tmesh = mesh;
 
  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    T tag = mesh.tag[eidx];
    
    for (int n = mesh.dsp[eidx]; n < mesh.dsp[eidx+1];n++)
    {
      T g = l2g[mesh.con[n]];
      T old_g = dof2vertex[g];
 
      const auto local_idx = g2l[g];
      const auto orig_idx  = g2l_orig[old_g];
 
      tmesh.xyz[local_idx*3+0] = Buffxyz_orig[orig_idx*3+0];
      tmesh.xyz[local_idx*3+1] = Buffxyz_orig[orig_idx*3+1];
      tmesh.xyz[local_idx*3+2] = Buffxyz_orig[orig_idx*3+2]; 
      dof2ptsData[g] = vertex2ptsdata[old_g];
 
      if (!std::isfinite(tmesh.xyz[local_idx*3+0]) || 
          !std::isfinite(tmesh.xyz[local_idx*3+1]) || 
          !std::isfinite(tmesh.xyz[local_idx*3+2]))
      {
        printf("NaN detected in coordinates of the EMI mesh after decoupling interfaces on EMI mesh while recording pts data on the EMI mesh, on Tag=%d\n",tag);
        EXIT(EXIT_FAILURE);
      }   
    }
  }
 
  mesh = tmesh;
}
 
template<class T, class S>
inline void insert_points_to_surface_mesh( const meshdata<T, S> & mesh_original, meshdata<T, S> & surfmesh, 
                                           const SF_nbr numbering,                             
                                           hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex,
                                           hashmap::unordered_set<int> & extra_tags,
                                           SF::vector<mesh_int_t> & elemTag_surface)
{
  MPI_Comm comm = mesh_original.comm;
 
  const SF::vector<mesh_int_t> & rnod_original = mesh_original.get_numbering(numbering);
  hashmap::unordered_map<T,T> g2l_orig;
  for(size_t i=0; i<rnod_original.size(); i++){
    g2l_orig[rnod_original[i]] = i;
  }
 
  const SF::vector<mesh_int_t> & rnod_surf = surfmesh.get_numbering(numbering);
 
  hashmap::unordered_map<T,T> g2l_surf;
  hashmap::unordered_map<T,T> l2g_surf;
  for(size_t i=0; i<rnod_surf.size(); i++){
    g2l_surf[rnod_surf[i]] = i;
    l2g_surf[i] = rnod_surf[i];
  }
 
  mesh_int_t gmax = global_max(rnod_surf, comm);
  surfmesh.g_numpts = gmax+1;
 
  SF::vector<mesh_real_t> Buffxyz = mesh_original.xyz;
 
  surfmesh.xyz.resize(surfmesh.l_numpts*3);
 
  elemTag_surface.resize(surfmesh.l_numelem);
  for(size_t eidx = 0; eidx < surfmesh.l_numelem; eidx++) {
    T tag = surfmesh.tag[eidx];
    elemTag_surface[eidx] = 1; // assigned 1 for for extracellular located on extracellular side
    if(extra_tags.find(tag) == extra_tags.end())
      elemTag_surface[eidx] = 2; // assigned 2 for for intracellular located on myocyte 
 
    for (int n = surfmesh.dsp[eidx]; n < surfmesh.dsp[eidx+1];n++)
    {
      T g = l2g_surf[surfmesh.con[n]];
      T old_g = dof2vertex[g];
 
      const auto local_idx = g2l_surf[g];
      const auto orig_idx  = g2l_orig[old_g];
 
      surfmesh.xyz[local_idx*3+0] = Buffxyz[orig_idx*3+0];
      surfmesh.xyz[local_idx*3+1] = Buffxyz[orig_idx*3+1];
      surfmesh.xyz[local_idx*3+2] = Buffxyz[orig_idx*3+2];
 
      if (!std::isfinite(surfmesh.xyz[local_idx*3+0]) || 
          !std::isfinite(surfmesh.xyz[local_idx*3+1]) || 
          !std::isfinite(surfmesh.xyz[local_idx*3+2]))
      {
        fprintf(stderr, "NaN detected in the surface mesh coordinates! on tag=%d \n", tag);
        exit(1);
      }
    }
  }
}
 
template<class T, class S>
inline void insert_points_to_surface_mesh( const meshdata<T, S> & mesh_original, meshdata<T, S> & surfmesh, 
                                           const SF_nbr numbering,                             
                                           hashmap::unordered_map<mesh_int_t,mesh_int_t> & dof2vertex)
{
  MPI_Comm comm = mesh_original.comm;
 
  const SF::vector<mesh_int_t> & rnod_original = mesh_original.get_numbering(numbering);
  hashmap::unordered_map<T,T> g2l_orig;
  for(size_t i=0; i<rnod_original.size(); i++){
    g2l_orig[rnod_original[i]] = i;
  }
 
  const SF::vector<mesh_int_t> & rnod_surf = surfmesh.get_numbering(numbering);
 
  hashmap::unordered_map<T,T> g2l_surf;
  hashmap::unordered_map<T,T> l2g_surf;
  for(size_t i=0; i<rnod_surf.size(); i++){
    g2l_surf[rnod_surf[i]] = i;
    l2g_surf[i] = rnod_surf[i];
  }
 
  mesh_int_t gmax = global_max(rnod_surf, comm);
  surfmesh.g_numpts = gmax+1;
 
  SF::vector<mesh_real_t> Buffxyz = mesh_original.xyz;
 
  surfmesh.xyz.resize(surfmesh.l_numpts*3);
 
  for(size_t eidx = 0; eidx < surfmesh.l_numelem; eidx++) {
    T tag = surfmesh.tag[eidx];
 
    for (int n = surfmesh.dsp[eidx]; n < surfmesh.dsp[eidx+1];n++)
    {
      T g = l2g_surf[surfmesh.con[n]];
      T old_g = dof2vertex[g];
 
      const auto local_idx = g2l_surf[g];
      const auto orig_idx  = g2l_orig[old_g];
 
      surfmesh.xyz[local_idx*3+0] = Buffxyz[orig_idx*3+0];
      surfmesh.xyz[local_idx*3+1] = Buffxyz[orig_idx*3+1];
      surfmesh.xyz[local_idx*3+2] = Buffxyz[orig_idx*3+2];
 
      if (!std::isfinite(surfmesh.xyz[local_idx*3+0]) || 
          !std::isfinite(surfmesh.xyz[local_idx*3+1]) || 
          !std::isfinite(surfmesh.xyz[local_idx*3+2]))
      {
        fprintf(stderr, "NaN detected in the surface mesh coordinates! on tag=%d \n", tag);
        exit(1);
      }
    }
  }
}
 
}
#endif
#endif