SF_numbering.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_NUMBERING_H
#define _SF_NUMBERING_H

#include <cstddef>
#include <mpi.h>

#include "hashmap.hpp"

#include "SF_container.h"
#include "SF_vector.h"

namespace SF {

template<class T, class S>
class numbering {
public:

  virtual void operator() (meshdata<T, S> & mesh) = 0;
};


template<class T, class S>
class submesh_numbering : public numbering<T, S>
{
public:
  inline void renumber_sorted_ascending(meshdata<T, S> & submesh, SF_nbr in_nbr, SF_nbr out_nbr)
  {
    MPI_Comm comm = submesh.comm;

    int size, rank;
    MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);

    // determine global min and max indices
    const vector<T> & rnod = submesh.get_numbering(in_nbr);
    T gmax = global_max(rnod, comm);
    T gmin = global_min(rnod, comm);

    // block size
    T bsize = (gmax - gmin) / size + 1;

    // distribute indices uniquely and linearly ascending (w.r.t. rank index) --------------------
    vector<T> dest(rnod.size()), sbuff(rnod);
    for(size_t i=0; i<dest.size(); i++)
      dest[i] = (rnod[i] - gmin) / bsize;
    binary_sort_copy(dest, sbuff);

    commgraph<size_t> grph;
    grph.configure(dest, comm);

    size_t rsize = sum(grph.rcnt);
    vector<T> rbuff(rsize);
    MPI_Exchange(grph, sbuff, rbuff, comm);

    // we need a data structure to store where we got the indices from
    vector<T> proc(rsize), acc_cnt(rsize, T(1));
    grph.source_ranks(proc);

    // compute local node set
    binary_sort_copy(rbuff, proc);
    unique_accumulate(rbuff, acc_cnt);
    rsize = rbuff.size();

    // communicate local sizes
    MPI_Allgather(&rsize, sizeof(size_t), MPI_BYTE, grph.rcnt.data(), sizeof(size_t), MPI_BYTE, comm);
    dsp_from_cnt(grph.rcnt, grph.rdsp);

    // now the new numbering can be derived from the data layout
    vector<T> new_nbr(rsize);
    interval(new_nbr, grph.rdsp[rank], grph.rdsp[rank+1]);

    // we send the old and new numberings back -------------------------------------------------
    sbuff.resize(proc.size());

    // expand the old numbering into the send buffer
    for(size_t i=0, idx=0; i<acc_cnt.size(); i++)
      for(T j=0; j<acc_cnt[i]; j++, idx++) sbuff[idx] = rbuff[i];
    dest.assign(proc.begin(), proc.end());
    binary_sort_copy(dest, sbuff); // permute for sending
    grph.configure(dest, comm);    // reconfigure comm graph
    rbuff.resize(sum(grph.rcnt));  // resize receive buffer
    MPI_Exchange(grph, sbuff, rbuff, comm);  // communicate

    vector<T> old_idx(rbuff);

    // expand the new numbering into the send buffer
    for(size_t i=0, idx=0; i<acc_cnt.size(); i++)
      for(T j=0; j<acc_cnt[i]; j++, idx++) sbuff[idx] = new_nbr[i];
    dest.assign(proc.begin(), proc.end());
    binary_sort_copy(dest, sbuff); // permute for sending
    MPI_Exchange(grph, sbuff, rbuff, comm);  // communicate

    vector<T> new_idx(rbuff);

    // add a new numbering to submesh -----------------------------------------
    // create mapping between original and new numbering
    index_mapping<T> map(old_idx, new_idx);
    vector<T> & snod = submesh.register_numbering(out_nbr);
    snod.resize(rnod.size());
    // use mapping for new numbering
    for(size_t i=0; i < rnod.size(); i++) snod[i] = map.forward_map(rnod[i]);
  }

  inline void operator() (meshdata<T, S> & submesh)
  {
    renumber_sorted_ascending(submesh, NBR_REF, NBR_SUBMESH);
    renumber_sorted_ascending(submesh, NBR_ELEM_REF, NBR_ELEM_SUBMESH);
  }

};

template<class T> inline
void reindex_cuthill_mckee(const vector<T> & n2n_cnt,
                           const vector<T> & n2n_dsp,
                           const vector<T> & n2n_con,
                           const bool reverse,
                           hashmap::unordered_map<T,T> & old2new)
{
  /*
   * We use the Cuthill-McKee algorithm to generate a permutation
   * of the original indices. This permutation is then renumbered
   * either forward or reverse.
   *
   * The permutation is generated as following:
   * Initially, an arbitrary node is selected. It is inserted in the permutation
   * and in the imaginary set R (represented by a boolean vector).
   * Then we loop over each node in perm, adding all connected nodes, that are not yet
   * in R, in ascending order to perm and R.
   *
   */
  size_t nnod = n2n_cnt.size();

  vector<bool>   inR(nnod, false);
  vector<T> perm(nnod);
  perm.resize(1); perm[0] = 0; inR[0] = true;
  T pidx=0;

  while(perm.size() < nnod)
  {
    T nidx;
    if (pidx < T(perm.size())) nidx = perm[pidx++];
    else {
      int i=0;
      while(inR[i] == true) i++;
      nidx = i;
      // std::cerr << "Warning: node " << nidx << " seems decoupled from mesh as it could not be reached by traversing the mesh edges!" << std::endl;
    }
    T start = n2n_dsp[nidx], stop = start + n2n_cnt[nidx];
    vector<T> adj(stop - start);

    T widx=0;
    for(T j = start; j<stop; j++) {
      T cidx = n2n_con[j];
      if( ! inR[cidx] ) {
        adj[widx++] = cidx;
        inR[cidx] = true;
      }
    }

    adj.resize(widx);
    binary_sort(adj);
    perm.append(adj.begin(), adj.end());
  }

  if(!reverse) {
    // Cuthill-McKee
    for(size_t i=0; i<perm.size(); i++) old2new[perm[i]] = i;
  } else {
    // Reverse Cuthill-McKee
    for(size_t i=0, j=perm.size()-1; i<perm.size(); i++, j--) old2new[perm[j]] = i;
  }
}


template<class T, class S>
class petsc_numbering : public numbering<T, S>
{
private:
  const overlapping_layout<T> & _pl;
public:
  petsc_numbering(const overlapping_layout<T> & pl) : _pl(pl)
  {}

  inline void operator() (meshdata<T, S> & mesh)
  {
    vector<T> & petsc_idx = mesh.register_numbering(NBR_PETSC);
    petsc_idx.assign(_pl.num_local_idx(), -1);

    int size, rank;
    MPI_Comm comm = mesh.comm;
    MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);


    const vector<T> & alg_dsp = _pl.algebraic_layout();
    const vector<T> & alg_nod = _pl.algebraic_nodes();
    const T alg_start = alg_dsp[rank], alg_end = alg_dsp[rank+1];

#if 1
    for(size_t i=0; i<alg_nod.size(); i++)
      petsc_idx[alg_nod[i]] = alg_start + i;
#else
    #if 0
    vector<bool> is_alg(mesh.l_numpts, false);
    for(const T & n : alg_nod) is_alg[n] = true;

    size_t loc_petsc_idx = 0;
    for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
      for(T j = mesh.dsp[eidx]; j < mesh.dsp[eidx+1]; j++) {
        const T c = mesh.con[j];
        if(is_alg[c] && petsc_idx[c] == -1) {
          petsc_idx[c] = alg_start + loc_petsc_idx;
          loc_petsc_idx++;
        }
      }
    }

    assert(loc_petsc_idx == alg_nod.size());
    #else
    meshdata<T,S> restr_mesh;
    vector<T> cnt(mesh.dsp.size() - 1, 0);
    hashmap::unordered_map<T,T> nodal2restr;
    restr_mesh.con.reserve(mesh.con.size());

    for(size_t i=0; i<alg_nod.size(); i++) {
      nodal2restr[alg_nod[i]] = i;
    }

    for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++)
    {
      T start = mesh.dsp[eidx], stop = mesh.dsp[eidx+1];
      for(T i = start; i<stop; i++) {
        if(nodal2restr.count(mesh.con[i])) {
          cnt[eidx]++;
          restr_mesh.con.push_back(nodal2restr[mesh.con[i]]);
        }
      }
    }

    dsp_from_cnt(cnt, restr_mesh.dsp);

    vector<T> n2n_cnt, n2n_dsp, n2n_con;
    hashmap::unordered_map<T,T> old2new;

    nodal_connectivity_graph(restr_mesh, n2n_cnt, n2n_con);
    dsp_from_cnt(n2n_cnt, n2n_dsp);

    reindex_cuthill_mckee(n2n_cnt, n2n_dsp, n2n_con, true, old2new);

    for(size_t i=0; i<alg_nod.size(); i++)
      petsc_idx[alg_nod[i]] = old2new[i] + alg_start;
    #endif
#endif
    // communicate the new numbering so that all nodes in the local DD domain
    // -- also the ones not in the algebraic range -- have a new index

    _pl.reduce(petsc_idx, "max");

    bool print_indexing = false;
    if(print_indexing) {
      for(int pid = 0; pid < size; pid++) {
        if(pid == rank) {
          for(size_t eidx = 0; eidx < 100; eidx++) {
            printf("rank %d elem %d: ", rank, int(eidx));
            for(T j = mesh.dsp[eidx]; j < mesh.dsp[eidx+1] - 1; j++) {
              const T c = mesh.con[j];
              printf("%d ", petsc_idx[c]);
            }
            printf("%d\n", petsc_idx[mesh.con[mesh.dsp[eidx+1] - 1]]);
          }
        }
        MPI_Barrier(PETSC_COMM_WORLD);
      }
    }
  }
};

}

#endif