SF_container.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
// ----------------------------------------------------------------------------
 
#ifndef _SF_CONTAINER_H
#define _SF_CONTAINER_H
 
#include <cstring>
#include <mpi.h>
#include <map>
#include <math.h>
#include <set>
 
#include "dense_mat.hpp"
#include "SF_globals.h"
#include "SF_vector.h"
#include "SF_linalg_utils.h"
#include "hashmap.hpp"
 
// This datatypes will be used in conjunction with mesh data. We might want to
// have increased / reduced precision when storing meshes, mesh related numberings, etc.
// Therefore, we will use (mesh_int_t,mesh_real_t) for data fields that scale with the mesh size,
// to be able to tweak the used datatype later. -Aurel
typedef int   mesh_int_t;
typedef float mesh_real_t;
 
namespace SF {
 
enum elem_t
{
  Tetra = 0,
  Hexa,
  Octa,
  Pyramid,
  Prism,
  Quad,
  Tri,
  Line
};
 
struct Point {
  double x;
  double y;
  double z;
};
 
template<class V>
inline Point arr_to_point(V* arr)
{
  return {arr[0], arr[1], arr[2]};
}
 
template<class V>
inline void point_to_arr(Point & p, V* arr)
{
  arr[0] = p.x, arr[1] = p.y, arr[2] = p.z;
}
 
inline Point cross(const Point & a, const Point & b)
{
  Point  c = {a.y*b.z - b.y*a.z, b.x*a.z - a.x*b.z, a.x*b.y - a.y*b.x};
  return c;
}
 
inline double inner_prod(const Point & a, const Point & b)
{
  return a.x*b.x + a.y*b.y + a.z*b.z;
}
 
inline void outer_prod(const Point & a, const Point & b, const double s, double* buff, const bool add = false)
{
  if(add == false)
    memset(buff, 0, 9*sizeof(double));
 
  buff[0] += a.x*b.x*s, buff[1] += a.x*b.y*s, buff[2] += a.x*b.z*s;
  buff[3] += a.y*b.x*s, buff[4] += a.y*b.y*s, buff[5] += a.y*b.z*s;
  buff[6] += a.z*b.x*s, buff[7] += a.z*b.y*s, buff[8] += a.z*b.z*s;
}
 
inline void outer_prod(const Point & a, const Point & b, double* buff)
{
  outer_prod(a, b, 1.0, buff);
}
 
inline double mag(const Point & vect)
{
  return sqrt(inner_prod(vect, vect));
}
 
inline Point normalize(const Point & vect)
{
  double norm = mag(vect);
  return {vect.x / norm, vect.y / norm, vect.z / norm};
}
 
inline Point operator-(const Point & a, const Point & b)
{
  return {a.x - b.x, a.y - b.y, a.z - b.z};
}
 
inline Point operator+(const Point & a, const Point & b)
{
  return {a.x + b.x, a.y + b.y, a.z + b.z};
}
 
inline Point operator*(const Point & a, const double & s)
{
  return {a.x * s, a.y * s, a.z * s};
}
 
inline Point operator*(const double & s, const Point & a)
{
  return {s * a.x, s * a.y, s * a.z};
}
 
inline Point operator*(const dmat<double>& A, const Point& b)
{
  return {A[0][0] * b.x + A[0][1] * b.y + A[0][2] * b.z,
          A[1][0] * b.x + A[1][1] * b.y + A[1][2] * b.z,
          A[2][0] * b.x + A[2][1] * b.y + A[2][2] * b.z};
}
 
inline Point project(const Point & a, const Point & b)
{
  return a * inner_prod(normalize(a), normalize(b));
}
 
inline Point orthogonalize(const Point & a, const Point & b)
{
  return a - project(a, b);
}
 
 
inline double distance(const Point & a, const Point & b)
{
  return mag(a - b);
}
 
inline elem_t getElemTypeID(char *eletype)
{
  elem_t ret;
 
  if ( !strcmp( eletype, "Tt" ) ) {
    ret = Tetra;
  } else if ( !strcmp( eletype, "Hx" ) ) {
    ret = Hexa;
  } else if ( !strcmp( eletype, "Oc" ) ) {
    ret = Octa;
  } else if ( !strcmp( eletype, "Py" ) ) {
    ret = Pyramid;
  } else if ( !strcmp( eletype, "Pr" ) ) {
    ret = Prism;
  } else if ( !strcmp( eletype, "Qd" ) ) {
    ret = Quad;
  } else if ( !strcmp( eletype, "Tr" ) ) {
    ret = Tri;
  } else {
    ret = Line;
  }
  return ret;
}
 
template<class S, class POINT> inline
void array_to_points(const vector<S> & arr, vector<POINT> & pts)
{
  pts.resize(arr.size() / 3);
  for(size_t i=0; i<pts.size(); i++)
    pts[i] = {arr[i*3+0], arr[i*3+1], arr[i*3+2]};
}
 
typedef enum {
  NBR_REF,           
  NBR_SUBMESH,       
  NBR_PETSC,         
  NBR_ELEM_REF,      
  NBR_ELEM_SUBMESH   
} SF_nbr;
 
// forward declaration of parallel layouts since they are used in mesh.
template<class T> class overlapping_layout;
template<class T> class non_overlapping_layout;
 
template<class T>
class index_mapping
{
  private:
  hashmap::unordered_map<T, T> _fwd;
  hashmap::unordered_map<T, T> _bwd;
 
  public:
  index_mapping()
  {}
 
  index_mapping(const vector<T> & a, const vector<T> & b)
  {
    this->assign(a, b);
  }
 
  inline void assign(const vector<T> & a, const vector<T> & b)
  {
    // a 1:1 mapping between old and new indices is required
    assert( a.size() == b.size() );
 
    _fwd.clear();
    _bwd.clear();
 
    for(size_t i=0; i<a.size(); i++)
    {
      T aidx = a[i], bidx = b[i];
      _fwd[aidx] = bidx;
      _bwd[bidx] = aidx;
    }
  }
 
  inline T forward_map(T idx) const
  {
    typename hashmap::unordered_map<T, T>::const_iterator it = _fwd.find(idx);
    if(it != _fwd.end())
      return it->second;
    else
      return T(-1);
  }
 
  inline T backward_map(T idx) const
  {
    typename hashmap::unordered_map<T, T>::const_iterator it = _bwd.find(idx);
    if(it != _bwd.end())
      return it->second;
    else
      return T(-1);
  }
 
 
  inline void forward_map(vector<T> & idx) const
  {
    size_t lsize = idx.size();
    vector<T> perm(lsize);
    interval(perm, 0, lsize);
 
    binary_sort_copy(idx, perm);
 
    size_t c = 0;
    typename hashmap::unordered_map<T, T>::const_iterator it;
 
    while(c < lsize)
    {
      T cc = idx[c];
      it = _fwd.find(cc);
 
      T val = it != _fwd.end() ? it->second : T(-1);
 
      while( (c < lsize) && (cc == idx[c]) )
      {
        idx[c] = val;
        c++;
      }
    }
 
    binary_sort_copy(perm, idx);
  }
 
  inline void backward_map(vector<T> & idx) const
  {
    size_t lsize = idx.size();
    vector<T> perm(lsize);
    interval(perm, 0, lsize);
 
    binary_sort_copy(idx, perm);
 
    size_t c = 0;
    typename hashmap::unordered_map<T, T>::const_iterator it;
 
    while(c < lsize)
    {
      T cc = idx[c];
      it = _bwd.find(cc);
 
      T val = it != _bwd.end() ? *it : T(-1);
 
      while( (c < lsize) && (cc == idx[c]) )
      {
        idx[c] = val;
        c++;
      }
    }
 
    binary_sort_copy(perm, idx);
  }
 
  size_t size() const
  {
    return _fwd.size();
  }
  bool in_a(const T idx) {
    return _fwd.count(idx) == 1;
  }
  bool in_b(const T idx) {
    return _bwd.count(idx) == 1;
  }
 
  const hashmap::unordered_map<T, T> & get_fwd_map() const
  {
    return _fwd;
  }
 
  const hashmap::unordered_map<T, T> & get_bwd_map() const
  {
    return _bwd;
  }
 
  hashmap::unordered_map<T, T> & get_fwd_map()
  {
    return _fwd;
  }
 
  hashmap::unordered_map<T, T> & get_bwd_map()
  {
    return _bwd;
  }
};
 
template<class T, class S>
class meshdata
{
public:
  size_t g_numelem; 
  size_t l_numelem; 
  size_t g_numpts;  
  size_t l_numpts;  
 
  MPI_Comm comm;
 
  std::string name;
 
  // element connectivity : the nodes forming the elements
  vector<T> con;
 
  // element data
  // one value per element
  vector<T> dsp;       
  vector<T> tag;       
  vector<elem_t> type; 
  // multiple values per element
  vector<S> fib;       
  vector<S> she;       
 
  hashmap::unordered_set<int> extr_tag;
 
  std::map< SF_nbr, vector<T> > nbr;   
  vector<S> xyz;                       
 
  overlapping_layout<T> pl;        
  non_overlapping_layout<T> epl;   
 
  meshdata(): g_numelem(0), l_numelem(0), g_numpts(0), l_numpts(0),
             comm(SF_COMM), name("unnamed")
  {}
 
  inline vector<T> & register_numbering(SF_nbr nbr_type)
  {
    if(nbr.count(nbr_type) == 0)
      nbr[nbr_type] = vector<T>();
 
    return nbr[nbr_type];
  }
 
  inline vector<T> & get_numbering(SF_nbr nbr_type)
  {
    typename std::map< SF_nbr, vector<T> >::iterator it = nbr.find(nbr_type);
    assert(it != nbr.end());
    return it->second;
  }
  inline const vector<T> & get_numbering(SF_nbr nbr_type) const
  {
    typename std::map< SF_nbr, vector<T> >::const_iterator it = nbr.find(nbr_type);
    assert(it != nbr.end());
    return it->second;
  }
 
  inline void localize(SF_nbr nbr_type)
  {
    // register numbering
    vector<T> & nod = this->register_numbering(nbr_type);
 
    hashmap::unordered_map<T,T> g2l;
    for(const T & c : this->con) g2l[c] = 0;
    g2l.sort();
 
    // compute numbering
    nod.resize(g2l.size());
 
    size_t widx = 0;
    for(auto it = g2l.begin(); it != g2l.end(); ++it, widx++) {
      nod[widx] = it->first;
      it->second = widx;
    }
 
    global_to_local(g2l, this->con, true);
    this->l_numpts = nod.size();
  }
 
  inline void globalize(SF_nbr nbr_type)
  {
    // register numbering
    vector<T> & nod = this->get_numbering(nbr_type);
    for(T & c : con) c = nod[c];
  }
 
 
  inline void generate_par_layout()
  {
    const vector<T> & rnod  = this->get_numbering(NBR_REF);
    const vector<T> & reidx = this->get_numbering(NBR_ELEM_REF);
 
    pl.assign(rnod, comm);
    g_numpts = pl.num_global_idx();
 
    epl.assign(reidx, comm);
  }
 
  inline void clear_data()
  {
    g_numelem = 0;
    g_numpts = 0;
    l_numelem = 0;
    l_numpts = 0;
 
    con.resize(0); con.reallocate();
    dsp.resize(0); dsp.reallocate();
    tag.resize(0); tag.reallocate();
    type.resize(0); type.reallocate();
 
    fib.resize(0); fib.reallocate();
    she.resize(0); she.reallocate();
    xyz.resize(0); xyz.reallocate();
 
    extr_tag.clear();
    nbr.clear();
  }
 
 
};
 
template<class T, class S>
inline void nodal_connectivity_graph(const meshdata<T, S> & mesh,
                              vector<T> & n2n_cnt,
                              vector<T> & n2n_con)
{
  vector<T> e2n_cnt, n2e_cnt, n2e_con;
  const vector<T> & e2n_con = mesh.con;
 
  cnt_from_dsp(mesh.dsp, e2n_cnt);
 
  transpose_connectivity(e2n_cnt, e2n_con, n2e_cnt, n2e_con);
  multiply_connectivities(n2e_cnt, n2e_con,
                          e2n_cnt, e2n_con,
                          n2n_cnt, n2n_con);
}
 
template<class T, class S>
int max_nodal_edgecount(const meshdata<T,S> & mesh)
{
  vector<T> n2n_cnt, n2n_con;
  nodal_connectivity_graph(mesh, n2n_cnt, n2n_con);
 
  int nmax = 0;
  for(size_t nidx=0; nidx < n2n_cnt.size(); nidx++)
    if(nmax < n2n_cnt[nidx]) nmax = n2n_cnt[nidx];
 
  MPI_Allreduce(MPI_IN_PLACE, &nmax, 1, MPI_INT, MPI_MAX, mesh.comm);
 
  return nmax;
}
 
template<class T, class S>
void get_alg_mask(const meshdata<T,S> & mesh, vector<bool> & alg_mask)
{
  const vector<T> & alg_nod = mesh.pl.algebraic_nodes();
 
  alg_mask.assign(mesh.l_numpts, false);
  for(const T & n : alg_nod) alg_mask[n] = true;
}
 
 
template<class T>
class commgraph
{
  public:
  vector<T> scnt; 
  vector<T> sdsp; 
  vector<T> rcnt; 
  vector<T> rdsp; 
 
  inline void resize(size_t size)
  {
    scnt.assign(size, T()), rcnt.assign(size, T()), sdsp.assign(size+1, T()), rdsp.assign(size+1, T());
  }
 
  template<class V>
  inline void scale(V fac)
  {
    for(size_t i=0; i<scnt.size(); i++)
      scnt[i] *= fac, rcnt[i] *= fac;
 
    for(size_t i=0; i<sdsp.size(); i++)
      sdsp[i] *= fac, rdsp[i] *= fac;
  }
 
  inline void transpose()
  {
    vector<T> tcnt(scnt), tdsp(sdsp);
    scnt = rcnt, sdsp = rdsp;
    rcnt = tcnt, rdsp = tdsp;
  }
 
  template<class V>
  inline void configure(const vector<V> & dest, MPI_Comm comm)
  {
    int size, rank;
    MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
    this->resize(size);
 
    count(dest, scnt);
    MPI_Alltoall(scnt.data(), sizeof(T), MPI_BYTE, rcnt.data(), sizeof(T), MPI_BYTE, comm);
    dsp_from_cnt(scnt, sdsp);
    dsp_from_cnt(rcnt, rdsp);
  }
  template<class V>
  inline void configure(const vector<V> & dest, const vector<T> & cnt, MPI_Comm comm)
  {
    int size, rank;
    MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);
 
    this->resize(size);
 
    scnt.zero();
    for(size_t i=0; i<dest.size(); i++) scnt[dest[i]] += cnt[i];
 
    MPI_Alltoall(scnt.data(), sizeof(T), MPI_BYTE, rcnt.data(), sizeof(T), MPI_BYTE, comm);
 
    dsp_from_cnt(scnt, sdsp);
    dsp_from_cnt(rcnt, rdsp);
  }
  template<class V>
  inline void source_ranks(vector<V> & source)
  {
    size_t size = rcnt.size();
    source.resize(sum(rcnt));
 
    for(size_t i=0, widx=0; i<size; i++)
      for(T j=0; j<rcnt[i]; j++, widx++)
        source[widx] = i;
  }
};
 
template<class T>
struct tuple {
  T v1;
  T v2;
};
 
template<class T, class S>
struct mixed_tuple
{
  T v1;
  S v2;
};
 
template<class T>
struct triple {
  T v1;
  T v2;
  T v3;
};
 
template<class T, class S, class V>
struct mixed_triple {
  T v1;
  S v2;
  V v3;
};
 
template<class T>
struct quadruple {
  T v1;
  T v2;
  T v3;
  T v4;
};
 
template<class T>
bool operator<(const struct tuple<T> & lhs, const struct tuple<T> & rhs)
{
  if(lhs.v1 != rhs.v1)
    return lhs.v1 < rhs.v1;
  else
    return lhs.v2 < rhs.v2;
}
template<class T, class S>
bool operator<(const struct mixed_tuple<T, S> & lhs, const struct mixed_tuple<T, S> & rhs)
{
  if(lhs.v1 != rhs.v1)
    return lhs.v1 < rhs.v1;
  else
    return lhs.v2 < rhs.v2;
}
 
template<class T>
bool operator<(const struct triple<T> & lhs, const struct triple<T> & rhs)
{
  if(lhs.v1 != rhs.v1)
    return lhs.v1 < rhs.v1;
  else if(lhs.v2 != rhs.v2)
    return lhs.v2 < rhs.v2;
  else
    return lhs.v3 < rhs.v3;
}
 
template<class T>
bool operator<(const struct quadruple<T> & lhs, const struct quadruple<T> & rhs)
{
  if(lhs.v1 != rhs.v1)
    return lhs.v1 < rhs.v1;
  else if(lhs.v2 != rhs.v2)
    return lhs.v2 < rhs.v2;
  else if(lhs.v3 != rhs.v3)
    return lhs.v3 < rhs.v3;
  else
    return lhs.v4 < rhs.v4;
}
 
// expand hashing for custom triple struct
 
}
 
namespace hashmap {
 
template<typename T>
struct hash_ops< SF::tuple<T> >
{
  static inline
  bool cmp(SF::tuple<T> a, SF::tuple<T> b) {
    return a.v1 == b.v1 && a.v2 == b.v2;
  }
 
  static inline
  hm_uint hash(SF::tuple<T> a) {
    hm_uint h = hashmap::mkhash_init;
    h = hashmap::mkhash(h, hashmap::hash_ops<T>::hash(a.v1));
    h = hashmap::mkhash(h, hashmap::hash_ops<T>::hash(a.v2));
    return h;
  }
};
 
template<typename T>
struct hash_ops< SF::triple<T> >
{
  static inline
  bool cmp(SF::triple<T> a, SF::triple<T> b) {
    return a.v1 == b.v1 && a.v2 == b.v2 && a.v3 == b.v3;
  }
 
  static inline
  hm_uint hash(SF::triple<T> a) {
    hm_uint h = mkhash_init;
    h = hashmap::mkhash(h, hashmap::hash_ops<T>::hash(a.v1));
    h = hashmap::mkhash(h, hashmap::hash_ops<T>::hash(a.v2));
    h = hashmap::mkhash(h, hashmap::hash_ops<T>::hash(a.v3));
    return h;
  }
};
 
template<typename T>
struct hash_ops< SF::quadruple<T> >
{
  static inline
  bool cmp(SF::quadruple<T> a, SF::quadruple<T> b) {
    return a.v1 == b.v1 && a.v2 == b.v2 && a.v3 == b.v3 && a.v4 == b.v4;
  }
 
  static inline
  hm_uint hash(SF::quadruple<T> a) {
    hm_uint h = mkhash_init;
    h = hashmap::mkhash(h, hashmap::hash_ops<T>::hash(a.v1));
    h = hashmap::mkhash(h, hashmap::hash_ops<T>::hash(a.v2));
    h = hashmap::mkhash(h, hashmap::hash_ops<T>::hash(a.v3));
    h = hashmap::mkhash(h, hashmap::hash_ops<T>::hash(a.v4));
    return h;
  }
};
}
 
template<class T>
struct tri_sele {
  T v1;
  T v2;
  T v3;
  T eidx;
};
 
template<class T>
struct quad_sele {
  T v1;
  T v2;
  T v3;
  T v4;
  T eidx;
};
 
template<class T>
void sort_triple(const T in1, const T in2, const T in3, T & out1, T & out2, T & out3)
{
    bool t12 = in1 < in2;
  bool t13 = in1 < in3;
  bool t23 = in2 < in3;
 
  if(t12) {
    //(123),(312),(132)
    if(t23) {
      //123
      out1=in1; out2=in2; out3=in3;
    }
    else {
      //(312),(132)
      if(t13) {
        //132
        out1=in1; out2=in3; out3=in2;
      }
      else {
        //312
        out1=in3; out2=in1; out3=in2;
      }
    }
  }
  else {
    //(213),(231),(321)
    if(t23) {
      //(213),(231)
      if(t13) {
        //213
        out1=in2; out2=in1; out3=in3;
      }
      else {
        //231
        out1=in2; out2=in3; out3=in1;
      }
    }
    else {
      //321
      out1=in3; out2=in2; out3=in1;
    }
  }
}
 
#endif