globalDOFNumbering.hh

#ifndef GLOBALDOFNUMBERING_HH
#define GLOBALDOFNUMBERING_HH
 
template<class X, class GFS>
class GlobalNumbering {
 
public:
 
  GlobalNumbering (const GFS& gfs)
   : my_rank(gfs.gridView().comm().rank()),
     ranks(gfs.gridView().comm().size()),
     claims(gfs, my_rank),
     global_index(gfs, -1.0)
  {
    using Dune::PDELab::Backend::native;
 
    // Construct Dirichlet constraints for processor boundaries; needed to not claim DOFs on boundary
    typedef typename GFS::template ConstraintsContainer<typename X::field_type>::Type CC;
    auto cc_proc_boundaries = CC();
    Dune::PDELab::NoDirichletConstraintsParameters pnbc;
    Dune::PDELab::constraints(pnbc,gfs,cc_proc_boundaries);
    Dune::PDELab::set_constrained_dofs(cc_proc_boundaries,-1.0,claims); // -1 processor boundaries
 
 
    // Communicate claims and retain maximum claimed rank per DOF
    Dune::PDELab::MaxDataHandle<GFS,X> maxdh(gfs,claims);
    gfs.gridView().communicate(maxdh,Dune::All_All_Interface,Dune::ForwardCommunication);
 
 
    // Iterate over DOFs, counting only those not claimed by higher ranks as local DOFs
    for (unsigned int i = 0; i < native(claims).N(); i++) {
 
      if (my_rank == native(claims)[i][0])
        local_dof_number++;
    }
 
    // Communicate local DOF numbers to allow consecutive global numbering
    std::vector<int> local_dof_numbers_(ranks);
    MPI_Allgather(&local_dof_number, 1, Dune::MPITraits<int>::getType(), local_dof_numbers_.data(), 1, Dune::MPITraits<int>::getType(), gfs.gridView().comm());
 
 
    // Compute offset within global numbering by adding DOF numbers of lower ranks
    for (int r = 0; r < ranks; r++) {
      global_dof_number += local_dof_numbers_[r];
      if (r < my_rank)
        global_dof_offset += local_dof_numbers_[r];
    }
 
    // Set vector to global DOF index, starting from offset and incrementing.
    // DOFs claimed by higher ranks are instead set to zero.
    int global_dof_counter = global_dof_offset;
    for (unsigned int i = 0; i < native(claims).N(); i++) {
 
      if (my_rank == native(claims)[i][0]) {
        native(global_index)[i][0] = global_dof_counter;
        global_dof_counter++;
      } else {
        native(global_index)[i][0] = 0.0;
      }
    }
 
    // Finally add together everything so that all DOFs claimed by others (currently at zero)
    // receive their global indices from neighbors
    Dune::PDELab::AddDataHandle<GFS,X> adddh(gfs,global_index);
    gfs.gridView().communicate(adddh,Dune::All_All_Interface,Dune::ForwardCommunication);
 
  }
 
  int claimedByRank(int localIndex) {
    using Dune::PDELab::Backend::native;
    return native(claims)[localIndex][0];
  }
 
  int numberOfLocalDOFs() {
    return local_dof_number;
  }
 
  int numberOfGlobalDOFs() {
    return global_dof_number;
  }
 
  int globalIndexOffset() {
    return global_dof_offset;
  }
 
  int globalIndex(int localIndex) {
    using Dune::PDELab::Backend::native;
    return native(global_index)[localIndex][0];
  }
 
  X globalIndexVector() {
    return global_index;
  }
 
  X claimsVector() {
    return claims;
  }
 
private:
  int my_rank;
  int ranks;
  X claims;
  X global_index;
  int local_dof_number = 0;
  int global_dof_number = 0;
  int global_dof_offset = 0;
 
};
 
#endif