SF_mesh_utils.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_MESH_UTILS_H
#define _SF_MESH_UTILS_H

#include <algorithm>
#include <iostream>
#include <string>

#include "asciiPlotter.hpp"

#include "SF_container.h"
#include "SF_network.h"
#include "SF_parallel_layout.h"
#include "SF_numbering.h"
#include "SF_sort.h"
#include "SF_mesh_io.h"
#include "SF_vector.h"

namespace SF {


template<class T, class S>
inline void permute_mesh(const meshdata<T, S> & inmesh, meshdata<T, S> & outmesh, const vector<T> & perm)
{
  assert(perm.size() == inmesh.l_numelem);

  outmesh.comm = inmesh.comm;

  outmesh.g_numelem = inmesh.g_numelem;
  outmesh.l_numelem = inmesh.l_numelem;

  // copy numberings
  outmesh.nbr = inmesh.nbr;

  const vector<T> & ref_eidx_in = inmesh.get_numbering(NBR_ELEM_REF);
  vector<T> & ref_eidx_out      = outmesh.get_numbering(NBR_ELEM_REF);

  size_t numelem = outmesh.l_numelem, numcon = inmesh.con.size();
  bool twoFib = inmesh.she.size() == (inmesh.l_numelem*3);

  outmesh.dsp .resize(numelem+1);
  outmesh.tag .resize(numelem);
  outmesh.type.resize(numelem);
  outmesh.fib .resize(numelem*3);
  if(twoFib)
    outmesh.she .resize(numelem*3);

  outmesh.con.resize(numcon);

  vector<T> cnt(numelem);
  T* elem = outmesh.con.data();

  for(size_t i=0; i<numelem; i++) {
    outmesh.tag[i]  = inmesh.tag[perm[i]];
    ref_eidx_out[i] = ref_eidx_in[perm[i]];
    outmesh.type[i] = inmesh.type[perm[i]];

    outmesh.fib[i*3+0] = inmesh.fib[perm[i]*3+0];
    outmesh.fib[i*3+1] = inmesh.fib[perm[i]*3+1];
    outmesh.fib[i*3+2] = inmesh.fib[perm[i]*3+2];

    if(twoFib) {
      outmesh.she[i*3+0] = inmesh.she[perm[i]*3+0];
      outmesh.she[i*3+1] = inmesh.she[perm[i]*3+1];
      outmesh.she[i*3+2] = inmesh.she[perm[i]*3+2];
    }

    int esize = inmesh.dsp[perm[i]+1] - inmesh.dsp[perm[i]];
    cnt[i] = esize;

    T estart = inmesh.dsp[perm[i]];
    for(int j=0; j<esize; j++) elem[j] = inmesh.con[estart+j];
    elem += esize;
  }

  dsp_from_cnt(cnt, outmesh.dsp);

  // check consistency
  assert(inmesh.dsp[inmesh.l_numelem] == outmesh.dsp[outmesh.l_numelem]);
}

template<class T, class S>
inline void redistribute_elements(meshdata<T, S> & mesh, meshdata<T, S> & sendbuff, vector<T> & part)
{
  MPI_Comm comm = mesh.comm;

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

  // first, we reorder the elements into contiguous blocks so we can
  // communicate them more easily
  vector<T> perm;
  interval(perm, 0, part.size());
  binary_sort_copy(part, perm);

  permute_mesh(mesh, sendbuff, perm);

  // now we set up the communication graphs (i.e. cnt and dsp for sending and receiving)
  commgraph<size_t> elem_grph, con_grph;
  // configure element exchange with the partitioning information.
  elem_grph.configure(part, comm);

  // configure connectivity exchange based on element exchange
  con_grph.resize(size);
  con_grph.sdsp[0] = 0;
  for(int i=0; i<size; i++) con_grph.sdsp[i+1] = sendbuff.dsp[elem_grph.sdsp[i+1]];
  cnt_from_dsp(con_grph.sdsp, con_grph.scnt);
  MPI_Alltoall(con_grph.scnt.data(), sizeof(size_t), MPI_BYTE, con_grph.rcnt.data(), sizeof(size_t), MPI_BYTE, comm);
  dsp_from_cnt(con_grph.rcnt, con_grph.rdsp);

  int twoFib = sendbuff.l_numelem > 0 ? sendbuff.she.size() == sendbuff.l_numelem*3 : 0;
  MPI_Allreduce(MPI_IN_PLACE, &twoFib, 1, MPI_INT, MPI_MAX, comm);

  // resize meshdata to fit received data
  vector<T> & ref_eidx       = mesh.get_numbering(NBR_ELEM_REF);
  vector<T> & ref_eidx_sbuff = sendbuff.get_numbering(NBR_ELEM_REF);

  size_t recv_size = sum(elem_grph.rcnt);
  mesh.l_numelem = recv_size;
  mesh.dsp .resize(recv_size+1);
  mesh.tag .resize(recv_size);
  ref_eidx.resize(recv_size);
  mesh.type.resize(recv_size);
  mesh.fib .resize(recv_size*3);
  if(twoFib)
    mesh.she .resize(recv_size*3);

  // we need auxiliary mesh counts
  vector<T> sendmesh_cnt(sendbuff.l_numelem), mesh_cnt(mesh.l_numelem);
  cnt_from_dsp(sendbuff.dsp, sendmesh_cnt);

  MPI_Exchange(elem_grph, sendmesh_cnt, mesh_cnt, comm);
  dsp_from_cnt(mesh_cnt, mesh.dsp);

  MPI_Exchange(elem_grph, sendbuff.tag,   mesh.tag,  comm);
  MPI_Exchange(elem_grph, ref_eidx_sbuff, ref_eidx,  comm);
  MPI_Exchange(elem_grph, sendbuff.type,  mesh.type, comm);

  elem_grph.scale(3);  // fiber data has three values per element
  MPI_Exchange(elem_grph, sendbuff.fib, mesh.fib, comm);
  if(twoFib)
    MPI_Exchange(elem_grph, sendbuff.she, mesh.she, comm);

  // the only thing left is to communicate connectivity
  {
    // we map the sendbuff connectivity to global reference indexing
    vector<T> & rnod = sendbuff.get_numbering(NBR_REF);
    for(size_t i=0; i<sendbuff.con.size(); i++) sendbuff.con[i] = rnod[sendbuff.con[i]];
  }
  recv_size = sum(con_grph.rcnt);
  mesh.con.resize(recv_size);
  MPI_Exchange(con_grph, sendbuff.con, mesh.con, comm);

  mesh.localize(NBR_REF);
}
template<class T, class S>
inline void redistribute_elements(meshdata<T, S> & mesh, vector<T> & part)
{
  meshdata<T, S> sendbuff;
  redistribute_elements(mesh, sendbuff, part);
}

template<class T, class S>
inline void redistribute_mesh(meshdata<T, S> & mesh, vector<T> & part)
{
  // the main idea is to compute the unique node set of the element block we send to
  // each process. Therefore we know which coordinates to communicate.

  meshdata<T, S> sendmesh;
  redistribute_elements(mesh, sendmesh, part);
  sendmesh.xyz.assign(mesh.xyz.begin(), mesh.xyz.end());
  // part and sendmesh have been sorted in redistribute_elements.

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

  // we use commgraphs to store the layouts of the elements, connectivities and nodes
  commgraph<T> elem_layout, con_layout, nod_layout;
  elem_layout.resize(size);
  con_layout.resize(size);
  nod_layout.resize(size);

  // compute layout of the elements
  count(part, elem_layout.scnt);
  dsp_from_cnt(elem_layout.scnt, elem_layout.sdsp);

  // compute layout of the connectivities using the element layout
  con_layout.scnt.zero();
  for(int pid=0; pid<size; pid++)
    for(T i=elem_layout.sdsp[pid]; i<elem_layout.sdsp[pid+1]; i++)
      con_layout.scnt[pid] += sendmesh.dsp[i+1] - sendmesh.dsp[i];
  dsp_from_cnt(con_layout.scnt, con_layout.sdsp);

  // finally compute the nodes we need to send to each rank and their layout
  vector<T> nod_sbuff(sendmesh.con.size());
  nod_layout.sdsp[0] = 0;  // initialize first entry of displacements

  for(int pid=0; pid<size; pid++)
  {
    T con_start = con_layout.sdsp[pid], con_end = con_layout.sdsp[pid+1];
    vector<T> nod(con_end - con_start);
    // copy connectivity block
    vec_assign(nod.data(), sendmesh.con.data() + con_start, nod.size());
    // get unique node set for block
    binary_sort(nod); unique_resize(nod);

    // add block size to layout
    nod_layout.scnt[pid]   = nod.size();
    nod_layout.sdsp[pid+1] = nod_layout.sdsp[pid] + nod_layout.scnt[pid];
    // copy node set to send buffer
    vec_assign(nod_sbuff.data() + nod_layout.sdsp[pid], nod.data(), nod.size());
  }
  nod_sbuff.resize(nod_layout.sdsp[nod_layout.sdsp.size()-1]);
  // finish up nod_layout
  MPI_Alltoall(nod_layout.scnt.data(), sizeof(T), MPI_BYTE, nod_layout.rcnt.data(), sizeof(T), MPI_BYTE, comm);
  dsp_from_cnt(nod_layout.rcnt, nod_layout.rdsp);

  vector<T> nod_lidx(nod_sbuff);
  global_to_local(sendmesh.get_numbering(NBR_REF), nod_lidx, false, true);

  // fill coordinates send buffer
  vector<S> xyz_sbuff(nod_sbuff.size() * 3);
  for(size_t i=0; i<nod_lidx.size(); i++)
  {
    T lidx = nod_lidx[i];
    xyz_sbuff[i*3+0] = sendmesh.xyz[lidx*3+0];
    xyz_sbuff[i*3+1] = sendmesh.xyz[lidx*3+1];
    xyz_sbuff[i*3+2] = sendmesh.xyz[lidx*3+2];
  }

  size_t rsize = sum(nod_layout.rcnt);
  vector<T> nod_rbuff(rsize);
  vector<S> xyz_rbuff(rsize*3);

  MPI_Exchange(nod_layout, nod_sbuff, nod_rbuff, comm);
  nod_layout.scale(3);
  MPI_Exchange(nod_layout, xyz_sbuff, xyz_rbuff, comm);

  // map received indices to local indexing
  global_to_local(mesh.get_numbering(NBR_REF), nod_rbuff, false, true);

  // find unique coordinates
  vector<T> acc_cnt(nod_rbuff.size(), 1), acc_dsp,
               acc_col(nod_rbuff.size());

  interval(acc_col, 0, nod_rbuff.size());
  binary_sort_copy(nod_rbuff, acc_col);
  unique_accumulate(nod_rbuff, acc_cnt);
  acc_dsp.resize(acc_cnt.size()+1);
  dsp_from_cnt(acc_cnt, acc_dsp);

  // generate new coordinates vector
  mesh.xyz.resize(acc_cnt.size()*3);
  for(size_t i=0; i<acc_cnt.size(); i++)
  {
    T pidx = acc_col[acc_dsp[i]];
    mesh.xyz[i*3+0] = xyz_rbuff[pidx*3+0];
    mesh.xyz[i*3+1] = xyz_rbuff[pidx*3+1];
    mesh.xyz[i*3+2] = xyz_rbuff[pidx*3+2];
  }
}


template<class T, class S>
inline void write_mesh_parallel(const meshdata<T, S> & mesh, bool binary, std::string basename)
{
  const MPI_Comm comm = mesh.comm;

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

  // we duplicate the mesh so that we can safely redistribute it
  meshdata<T, S> wmesh = mesh;
  vector<T> & ref_eidx = wmesh.get_numbering(NBR_ELEM_REF);

  T elem_gmin = global_min(ref_eidx, comm);
  T elem_gmax = global_max(ref_eidx, comm);

  // block size
  T elem_bsize = (elem_gmax - elem_gmin) / size + 1;

  // redistribute elements linearly ascending after their reference element index --------------------
  vector<T> dest(wmesh.l_numelem);
  for(size_t i=0; i<dest.size(); i++)
    dest[i] = (ref_eidx[i] - elem_gmin) / elem_bsize;

  // redistribute elements of write mesh.
  // WARNING: vertex coordinates are not redistributed. therefore the
  // vertices will be redistributed based on the original mesh.
  redistribute_elements(wmesh, dest);

  {
    // temporary datastructs
    vector<T> perm, eidx;
    meshdata<T, S> tmesh = wmesh;

    // generate permutation
    eidx = wmesh.get_numbering(NBR_ELEM_REF);
    interval(perm, 0, eidx.size());
    binary_sort_copy(eidx, perm);

    // permute local mesh to locally ascending element indices
    permute_mesh(tmesh, wmesh, perm);
  }

  vector<T> nbr_orig; // the global vertex numbering of the original mesh
  if(mesh.nbr.count(NBR_SUBMESH)) {
    // mesh is a submesh. generate submesh numbering for the redistributed
    // mesh and map connectivity to it

    // copy original numbering for later
    nbr_orig = mesh.get_numbering(NBR_SUBMESH);

    // generate new numbering on write mesh
    submesh_numbering<T, S> submesh_nbr;
    submesh_nbr(wmesh);
    const vector<T> & nbr = wmesh.get_numbering(NBR_SUBMESH);
    for(size_t i=0; i<wmesh.con.size(); i++) wmesh.con[i] = nbr[wmesh.con[i]];
  }
  else {
    // mesh is a reference mesh. just map connectivity.

    // copy original numbering for later
    nbr_orig = mesh.get_numbering(NBR_REF);

    const vector<T> & nbr = wmesh.get_numbering(NBR_REF);
    for(size_t i=0; i<wmesh.con.size(); i++) wmesh.con[i] = nbr[wmesh.con[i]];
  }

  // write .elem and .lon files
  write_elements(wmesh, binary, basename);

  // now we need to redistribute and write the vertex coordinates
  const vector<T> & alg_nod = mesh.pl.algebraic_nodes();
  vector<T> xyz_cnt(alg_nod.size(), 3), xyz_idx(alg_nod.size()),  srt_cnt, srt_idx;
  vector<S> xyz(alg_nod.size()*3), srt_xyz;

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

    xyz[i*3+0] = mesh.xyz[loc*3+0];
    xyz[i*3+1] = mesh.xyz[loc*3+1];
    xyz[i*3+2] = mesh.xyz[loc*3+2];

    xyz_idx[i] = nbr_orig[loc];
  }

  sort_parallel(comm, xyz_idx, xyz_cnt, xyz, srt_idx, srt_cnt, srt_xyz);

  // write header with root
  FILE* pts_fd;
  std::string pts_file = binary ? basename + ".bpts" : basename + ".pts";

  if(rank == 0) {
    pts_fd = fopen(pts_file.c_str(), "w");
    if(!pts_fd) {
      fprintf(stderr, "Error: could not open file: %s. Aborting!", pts_file.c_str());
      exit(1);
    }
    write_pts_header(pts_fd, binary, mesh.g_numpts);
    fclose(pts_fd);
  }

  // write vertices sequentially rank after rank
  for(int pid=0; pid < size; pid++) {
    if(pid == rank) {
      pts_fd = fopen(pts_file.c_str(), "a");
      if(!pts_fd) {
        fprintf(stderr, "Error: could not open file: %s. Aborting!", pts_file.c_str());
        exit(1);
      }
      write_pts_block(pts_fd, binary, srt_xyz);
      fclose(pts_fd);
    }
    MPI_Barrier(comm);
  }
}


template<class T, class S>
inline void gather_mesh(const meshdata<T, S> & locmesh, meshdata<T, S> & globmesh)
{
  MPI_Comm comm = locmesh.comm;

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

  // gather mesh data using the redistribute_mesh function
  vector<T> part(locmesh.l_numelem, 0); // all elements to rank 0
  globmesh = locmesh;
  redistribute_mesh(globmesh, part);
}

template<class T, class S>
inline void print_DD_info(const meshdata<T, S> & mesh)
{
  MPI_Comm comm = mesh.comm;

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

  vector<size_t> npoint(size), nelem(size), ninterf(size), nidx(size);

  size_t intf_size = mesh.pl.interface().size(),
         idx_size  = mesh.pl.num_algebraic_idx();

  MPI_Gather(&mesh.l_numpts,  sizeof(size_t), MPI_BYTE, npoint.data(),  sizeof(size_t), MPI_BYTE, 0, comm);
  MPI_Gather(&mesh.l_numelem, sizeof(size_t), MPI_BYTE, nelem.data(),   sizeof(size_t), MPI_BYTE, 0, comm);
  MPI_Gather(&intf_size,      sizeof(size_t), MPI_BYTE, ninterf.data(), sizeof(size_t), MPI_BYTE, 0, comm);
  MPI_Gather(&idx_size,       sizeof(size_t), MPI_BYTE, nidx.data(),    sizeof(size_t), MPI_BYTE, 0, comm);

  vector<int> mult(mesh.l_numpts, 1);
  mesh.pl.reduce(mult, "sum");

  int hist_size = 16;
  vector<int> mult_hist(hist_size, 0), global_mult_hist(hist_size, 0);
  for(auto m : mult) mult_hist[m]++;

  MPI_Reduce(mult_hist.data(), global_mult_hist.data(), hist_size, MPI_INT, MPI_SUM, 0, comm);

  if(!rank) {
    printf("===== Parallel mesh statistics =====\n");

    printf("#pid\t#nodes\t#elems\t#interf\t#alg\n");
    for(int pid = 0; pid < size; pid++)
      printf("%d\t%ld\t%ld\t%ld\t%ld\n", pid, (long int)npoint[pid], (long int)nelem[pid],
             (long int)ninterf[pid], (long int)nidx[pid]);
    printf("\n");

    std::cout << "Multiplicities :" << std::endl;
    for(int i = 2; i < hist_size && global_mult_hist[i] > 0; i++)
      std::cout << i << ": " << global_mult_hist[i] << std::endl;
  }
}


template<class T, class S>
inline void extract_mesh(const vector<bool> & keep,
                         const meshdata<T, S> & mesh,
                         meshdata<T, S> & submesh)
{
  assert(keep.size() == mesh.l_numelem);

  size_t num_extr_elem = 0, num_extr_entr = 0;
  bool twoFib = ( mesh.she.size() == (mesh.l_numelem*3) );

  submesh.l_numpts = submesh.g_numpts = 0;
  submesh.comm = mesh.comm;

  const vector<T> & mesh_ref_eidx = mesh.get_numbering(NBR_ELEM_REF);
  vector<T> & sub_ref_eidx        = submesh.register_numbering(NBR_ELEM_REF);

  for(size_t i=0; i<mesh.l_numelem; i++) {
    if(keep[i]) {
      num_extr_elem++;
      num_extr_entr += mesh.dsp[i+1] - mesh.dsp[i];
    }
  }

  submesh.l_numelem = num_extr_elem;

  vector<T> cnt(num_extr_elem);
  submesh.dsp.resize(num_extr_elem+1);
  submesh.tag.resize(num_extr_elem);
  sub_ref_eidx.resize(num_extr_elem);
  submesh.type.resize(num_extr_elem);

  submesh.fib.resize(num_extr_elem*3);
  if(twoFib)
    submesh.she.resize(num_extr_elem*3);

  submesh.con.resize(num_extr_entr);
  // reference numbering
  const vector<T> & rnod = mesh.get_numbering(NBR_REF);

  for(size_t ridx_ele=0, ridx_con=0, widx_ele=0, widx_con=0; ridx_ele<mesh.l_numelem; ridx_ele++)
  {
    ridx_con = mesh.dsp[ridx_ele];

    if(keep[ridx_ele]) {
      // element data
      cnt[widx_ele] = mesh.dsp[ridx_ele + 1] - mesh.dsp[ridx_ele];
      submesh.tag [widx_ele] = mesh.tag     [ridx_ele];
      sub_ref_eidx[widx_ele] = mesh_ref_eidx[ridx_ele];
      submesh.type[widx_ele] = mesh.type    [ridx_ele];

      // fibers
      submesh.fib[widx_ele*3+0] = mesh.fib[ridx_ele*3+0];
      submesh.fib[widx_ele*3+1] = mesh.fib[ridx_ele*3+1];
      submesh.fib[widx_ele*3+2] = mesh.fib[ridx_ele*3+2];
      if(twoFib) {
        submesh.she[widx_ele*3+0] = mesh.she[ridx_ele*3+0];
        submesh.she[widx_ele*3+1] = mesh.she[ridx_ele*3+1];
        submesh.she[widx_ele*3+2] = mesh.she[ridx_ele*3+2];
      }

      // connectivity
      for(int j=0; j<cnt[widx_ele]; j++)
        submesh.con[widx_con++] = rnod[mesh.con[ridx_con++]];

      widx_ele++;
    }
  }
  dsp_from_cnt(cnt, submesh.dsp);

  unsigned long int gnumele = submesh.l_numelem;
  MPI_Allreduce(MPI_IN_PLACE, &gnumele, 1, MPI_UNSIGNED_LONG, MPI_SUM, submesh.comm);
  submesh.g_numelem = gnumele;

  submesh.localize(NBR_REF);
}


template<class T, class S>
inline void rebalance_mesh(meshdata<T, S> & mesh)
{
  MPI_Comm comm = mesh.comm;

  // we check if some ranks have 0 local elems
  short errflag = 0;
  if(mesh.l_numelem == 0) errflag = 1;
  MPI_Allreduce(MPI_IN_PLACE, &errflag, 1, MPI_SHORT, MPI_SUM, comm);

  // we redistribute the elements if necessary
  if(errflag) {
    int size, rank;
    MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);

    vector<T> part(mesh.l_numelem);
    vector<size_t> elem_counts(size), layout;
    MPI_Allgather(&mesh.l_numelem, sizeof(size_t), MPI_BYTE, elem_counts.data(),
                  sizeof(size_t), MPI_BYTE, comm);

    dsp_from_cnt(elem_counts, layout);

    for(size_t i=0; i<mesh.l_numelem; i++)
      part[i] = (layout[rank] + i) % size;

    redistribute_elements(mesh, part);
  }
}


template<class T, class S>
inline void extract_myocardium(const meshdata<T, S> & mesh, meshdata<T, S> & submesh)
{
  MPI_Comm comm = mesh.comm;
  vector<bool> keep(mesh.l_numelem);

  // fill keep based on empty fiber definitions
  for(size_t i=0; i<mesh.l_numelem; i++) {
    S l1 = mesh.fib[i*3+0];
    S l2 = mesh.fib[i*3+1];
    S l3 = mesh.fib[i*3+2];

    if( l1*l1 + l2*l2 + l3*l3 )
      keep[i] = true;
    else
      keep[i] = false;
  }

  extract_mesh(keep, mesh, submesh);
  rebalance_mesh(submesh);

  // get unique tag set
  submesh.extr_tag.insert(submesh.tag.begin(), submesh.tag.end());
  {
    // we have to collect the global list of tags of this submesh
    vector<T> tags, glob_tags;
    tags.assign(submesh.extr_tag.begin(), submesh.extr_tag.end()); // local list
    make_global(tags, glob_tags, mesh.comm);                       // now globalized

    submesh.extr_tag.insert(glob_tags.begin(), glob_tags.end());
  }
}


template<class T, class S>
inline void extract_tagbased(const meshdata<T, S> & mesh, meshdata<T, S> & submesh)
{
  if(submesh.extr_tag.size())
  {
    vector<bool> keep(mesh.l_numelem, false);

    // fill keep based on empty fiber definitions
    for(size_t i=0; i<mesh.l_numelem; i++) {
      if( submesh.extr_tag.count(mesh.tag[i]) )
        keep[i] = true;
    }

    extract_mesh(keep, mesh, submesh);
  }
  else {
    // if the tags are not provided, we assume that the whole mesh is copied into the
    // submesh
    std::string oldname = submesh.name;
    // copy reference mesh
    submesh = mesh;
    // revert name to old name
    submesh.name = oldname;
    // compute tags
    submesh.extr_tag.insert(submesh.tag.begin(), submesh.tag.end());
    {
      // we have to collect the global list of tags of this submesh
      vector<T> tags, glob_tags;
      tags.assign(submesh.extr_tag.begin(), submesh.extr_tag.end()); // local list
      make_global(tags, glob_tags, mesh.comm);                       // now globalized

      submesh.extr_tag.insert(glob_tags.begin(), glob_tags.end());
    }
  }

  rebalance_mesh(submesh);
}

template<class T, class S>
inline void print_mesh_graph(meshdata<T,S> & mesh)
{
  int rank, size;
  MPI_Comm_size(mesh.comm, &size); MPI_Comm_rank(mesh.comm, &rank);

  vector<T> n2n_cnt, n2n_con;
  nodal_connectivity_graph(mesh, n2n_cnt, n2n_con);
  matrixgraph_plotter<vector<T> > plotter(30, 60);

  for(int pid=0; pid<size; pid++)
  {
    if(pid == rank) {
      std::cout << "\n\n Rank " << rank << ": \n" << std::endl;
      plotter.print(n2n_cnt, n2n_con, '*');
    }
    MPI_Barrier(mesh.comm);
  }
}

template<class T, class S> inline
void inter_domain_mapping(const meshdata<T,S> & mesh_a, const meshdata<T,S> & mesh_b,
                          const SF_nbr snbr, index_mapping<T> & a_to_b)
{
  assert(mesh_a.comm == mesh_b.comm);
  MPI_Comm comm = mesh_a.comm;

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

  // reference and submesh numberings of both meshes
  const vector<T> & mesh_a_rnbr = mesh_a.get_numbering(NBR_REF);
  const vector<T> & mesh_a_snbr = mesh_a.get_numbering(snbr);
  const vector<T> & mesh_b_rnbr = mesh_b.get_numbering(NBR_REF);
  const vector<T> & mesh_b_snbr = mesh_b.get_numbering(snbr);
  // send- and receive-buffer for the reference and submesh numberings
  vector<T> mesh_a_rnbr_sbuff, mesh_a_rnbr_rbuff, mesh_a_snbr_sbuff, mesh_a_snbr_rbuff;
  vector<T> mesh_b_rnbr_sbuff, mesh_b_rnbr_rbuff, mesh_b_snbr_sbuff, mesh_b_snbr_rbuff;

  // communicate reference- and submesh numbering of mesh a ----------------------------
  vector<T> dest(mesh_a_rnbr.size()), perm_a, perm_b;
  for(size_t i=0; i<dest.size(); i++) dest[i] = mesh_a_rnbr[i] % size;

  interval(perm_a, 0, mesh_a_rnbr.size());
  binary_sort_copy(dest, perm_a);

  mesh_a_rnbr_sbuff.resize(dest.size()); mesh_a_snbr_sbuff.resize(dest.size());
  for(size_t i=0; i<dest.size(); i++) {
    mesh_a_rnbr_sbuff[i] = mesh_a_rnbr[perm_a[i]];
    //mesh_a_snbr_sbuff[i] = mesh_a_snbr[perm_a[i]];
  }

  commgraph<size_t> grph_a;
  grph_a.configure(dest, comm);
  size_t rcv_size = sum(grph_a.rcnt);

  mesh_a_rnbr_rbuff.resize(rcv_size);
  mesh_a_snbr_rbuff.resize(rcv_size);

  MPI_Exchange(grph_a, mesh_a_rnbr_sbuff, mesh_a_rnbr_rbuff, comm);
  //MPI_Exchange(grph_a, mesh_a_snbr_sbuff, mesh_a_snbr_rbuff, comm);

  // communicate reference- and submesh numbering of mesh b ----------------------------
  dest.resize(mesh_b_rnbr.size());
  for(size_t i=0; i<dest.size(); i++) dest[i] = mesh_b_rnbr[i] % size;

  interval(perm_b, 0, mesh_b_rnbr.size());
  binary_sort_copy(dest, perm_b);

  mesh_b_rnbr_sbuff.resize(dest.size()); mesh_b_snbr_sbuff.resize(dest.size());
  for(size_t i=0; i<dest.size(); i++) {
    mesh_b_rnbr_sbuff[i] = mesh_b_rnbr[perm_b[i]];
    mesh_b_snbr_sbuff[i] = mesh_b_snbr[perm_b[i]];
  }

  commgraph<size_t> grph_b;
  grph_b.configure(dest, comm);
  rcv_size = sum(grph_b.rcnt);

  mesh_b_rnbr_rbuff.resize(rcv_size);
  mesh_b_snbr_rbuff.resize(rcv_size);

  MPI_Exchange(grph_b, mesh_b_rnbr_sbuff, mesh_b_rnbr_rbuff, comm);
  MPI_Exchange(grph_b, mesh_b_snbr_sbuff, mesh_b_snbr_rbuff, comm);

  // generate mapping ------------------------------------------------------------------
  // We do this by first setting up a map between reference index and submesh index for
  // mesh b. This allows us to check fast if a reference index is in b and do the mapping.
  hashmap::unordered_map<T,T> ref_to_sub_b;
  for(size_t i=0; i<mesh_b_rnbr_rbuff.size(); i++) {
    T ref_idx = mesh_b_rnbr_rbuff[i];
    T sub_idx = mesh_b_snbr_rbuff[i];

    if(ref_to_sub_b.count(ref_idx) && ref_to_sub_b[ref_idx] != sub_idx)
      fprintf(stderr, "inter_domain_mapping error: Missmatching multiple mappings: %d : %d \n",
              ref_to_sub_b[ref_idx], sub_idx);
    ref_to_sub_b[ref_idx] = sub_idx;
  }

  // Now we either map the reference indices of mesh a to the submesh index of b, or, if
  // the reference index is not in b, set it to -1
  for(size_t i=0; i<mesh_a_rnbr_rbuff.size(); i++) {
    auto it = ref_to_sub_b.find(mesh_a_rnbr_rbuff[i]);
    if(it != ref_to_sub_b.end())
      mesh_a_snbr_rbuff[i] = it->second;
    else
      mesh_a_snbr_rbuff[i] = -1;
  }

  // now we simply communicate back by transposing grph_a
  grph_a.transpose();
  MPI_Exchange(grph_a, mesh_a_snbr_rbuff, mesh_a_snbr_sbuff, comm);

  size_t num_mapped = 0;
  for(size_t i=0; i<mesh_a_snbr_sbuff.size(); i++)
    if(mesh_a_snbr_sbuff[i] > -1) num_mapped++;

  // the submesh indices that could be mapped
  vector<T> snbr_a(num_mapped), snbr_b(num_mapped);
  binary_sort_copy(perm_a, mesh_a_snbr_sbuff);

  for(size_t i=0, idx=0; i<mesh_a_snbr_sbuff.size(); i++) {
    if(mesh_a_snbr_sbuff[i] > -1) {
      snbr_a[idx] = mesh_a_snbr[i];
      snbr_b[idx] = mesh_a_snbr_sbuff[i];
      idx++;
    }
  }

  // now we can set up the a-to-b mapping
  a_to_b.assign(snbr_a, snbr_b);
}

template<class T> inline
void insert_surf_tri(T n1, T n2, T n3, size_t eidx,
                     triple<T> & surf,
                     tri_sele<T> & sele,
                     hashmap::unordered_map<triple<T>, tri_sele<T> > & surfmap)
{
  sele.v1 = n1, sele.v2 = n2, sele.v3 = n3; sele.eidx = eidx;
  sort_triple(n1, n2, n3, surf.v1, surf.v2, surf.v3);

  auto it = surfmap.find(surf);
  if(it != surfmap.end()) surfmap.erase(it);
  else surfmap[surf] = sele;
}

template<class T> inline
void insert_surf_quad(T n1, T n2, T n3, T n4, size_t eidx,
                      vector<T> & buff,
                      quadruple<T> & surf,
                      quad_sele<T> & sele,
                      hashmap::unordered_map<quadruple<T>, quad_sele<T> > & surfmap)
{
  buff[0] = n1, buff[1] = n2, buff[2] = n3, buff[3] = n4;
  binary_sort(buff);

  sele.v1 = n1, sele.v2 = n2, sele.v3 = n3, sele.v4 = n4, sele.eidx = eidx;
  surf.v1 = buff[0], surf.v2 = buff[1], surf.v3 = buff[2], surf.v4 = buff[3];

  auto it = surfmap.find(surf);
  if(it != surfmap.end()) surfmap.erase(it);
  else surfmap[surf] = sele;
}

template<class T> inline
void insert_surf_tet(const T* nod,
                     const size_t eidx,
                     hashmap::unordered_map<triple<T>, tri_sele<T> > & surfmap)
{
  // surfaces are (2,3,1) , (1,4,2) , (2,4,3) , (1,3,4)
  struct triple<T>   surf;
  struct tri_sele<T> sele;

  T n1 = nod[0], n2 = nod[1], n3 = nod[2], n4 = nod[3];

  insert_surf_tri(n2, n3, n1, eidx, surf, sele, surfmap);
  insert_surf_tri(n1, n4, n2, eidx, surf, sele, surfmap);
  insert_surf_tri(n2, n4, n3, eidx, surf, sele, surfmap);
  insert_surf_tri(n1, n3, n4, eidx, surf, sele, surfmap);
}

template<class T> inline
void insert_surf_pyr(const T* nod, const size_t eidx, vector<T> & buff,
                     hashmap::unordered_map<triple<T>, tri_sele<T> > & surfmap,
                     hashmap::unordered_map<quadruple<T>, quad_sele<T> > & qsurfmap)
{
  // surfaces are (1,5,2) , (2,5,3) , (3,5,4) , (4,5,1) , (1,2,3,4)
  triple<T>   surf; quadruple<T> qsurf;
  tri_sele<T> sele; quad_sele<T> qsele;

  T n1 = nod[0], n2 = nod[1], n3 = nod[2], n4 = nod[3], n5 = nod[4];

  insert_surf_tri(n1, n5, n2, eidx, surf, sele, surfmap);
  insert_surf_tri(n2, n5, n3, eidx, surf, sele, surfmap);
  insert_surf_tri(n3, n5, n4, eidx, surf, sele, surfmap);
  insert_surf_tri(n4, n5, n1, eidx, surf, sele, surfmap);
  insert_surf_quad(n1, n2, n3, n4, eidx, buff, qsurf, qsele, qsurfmap);
}

template<class T> inline
void insert_surf_pri(const T* nod, const size_t eidx, vector<T> & buff,
                     hashmap::unordered_map<triple<T>, tri_sele<T> > & surfmap,
                     hashmap::unordered_map<quadruple<T>, quad_sele<T> > & qsurfmap)
{
  struct triple<T>    surf;
  struct quadruple<T> qsurf;
  struct tri_sele<T>  sele;
  struct quad_sele<T> qsele;

  T n1 = nod[0], n2 = nod[1], n3 = nod[2], n4 = nod[3], n5 = nod[4], n6 = nod[5];

  // surfaces are (1,2,3) , (4,5,6) , (1,2,6,4) , (2,3,5,6) , (3,1,4,5)
  insert_surf_tri(n1, n2, n3, eidx, surf, sele, surfmap);
  insert_surf_tri(n4, n5, n6, eidx, surf, sele, surfmap);
  insert_surf_quad(n1, n2, n6, n4, eidx, buff, qsurf, qsele, qsurfmap);
  insert_surf_quad(n2, n3, n5, n6, eidx, buff, qsurf, qsele, qsurfmap);
  insert_surf_quad(n3, n1, n4, n5, eidx, buff, qsurf, qsele, qsurfmap);
}

template<class T> inline
void insert_surf_hex(const T* nod, const size_t eidx, vector<T> & buff,
                     hashmap::unordered_map<quadruple<T>, quad_sele<T> > & surfmap)
{
  // surfaces are (1,2,3,4) , (3,2,8,7) , (4,3,7,6) , (1,4,6,5) , (2,1,5,8), (6,7,8,5)
  struct quadruple<T> qsurf;
  struct quad_sele<T> qsele;

  T n1 = nod[0], n2 = nod[1], n3 = nod[2], n4 = nod[3],
    n5 = nod[4], n6 = nod[5], n7 = nod[6], n8 = nod[7];

  insert_surf_quad(n1, n2, n3, n4, eidx, buff, qsurf, qsele, surfmap);
  insert_surf_quad(n3, n2, n8, n7, eidx, buff, qsurf, qsele, surfmap);
  insert_surf_quad(n4, n3, n7, n6, eidx, buff, qsurf, qsele, surfmap);
  insert_surf_quad(n1, n4, n6, n5, eidx, buff, qsurf, qsele, surfmap);
  insert_surf_quad(n2, n1, n5, n8, eidx, buff, qsurf, qsele, surfmap);
  insert_surf_quad(n6, n7, n8, n5, eidx, buff, qsurf, qsele, surfmap);
}

template<class T, class S> inline
void search_for_surface(const meshdata<T,S> & mesh, const SF_nbr numbering,
                        const hashmap::unordered_map<triple<T>, tri_sele<T> >     & search_tri,
                        const hashmap::unordered_map<quadruple<T>, quad_sele<T> > & search_quad,
                        hashmap::unordered_map<triple<T>, tri_sele<T> >     & found_tri,
                        hashmap::unordered_map<quadruple<T>, quad_sele<T> > & found_quad)
{
  const T* con = mesh.con.data();
  const T* nbr = mesh.get_numbering(numbering).data();

  hashmap::unordered_map<triple<T>, tri_sele<T> >     tri_surf;
  hashmap::unordered_map<quadruple<T>, quad_sele<T> > quad_surf;

  vector<T> nodvec(mesh.con.size()); T* nod = nodvec.data();
  for(size_t i=0; i<nodvec.size(); i++)
   nod[i] = nbr[con[i]];

  const vector<T> & ref_eidx = mesh.get_numbering(NBR_ELEM_REF);
  vector<T> qbuff(4); // buffer if quad nodes need to be sorted

  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    tri_surf.clear();
    quad_surf.clear();

    switch(mesh.type[eidx]) {
      case Tri: {
        triple<T> tri;
        tri_sele<T> trie;
        insert_surf_tri(nod[0], nod[1], nod[2], ref_eidx[eidx], tri, trie, tri_surf);
        break;
      }

      case Quad: {
        quadruple<T> qd;
        quad_sele<T> qde;
        insert_surf_quad(nod[0], nod[1], nod[2], nod[3],
                         ref_eidx[eidx], qbuff, qd, qde, quad_surf);
        break;
      }

      case Tetra:
        insert_surf_tet(nod, ref_eidx[eidx], tri_surf);
        break;

      case Pyramid:
        insert_surf_pyr(nod, ref_eidx[eidx], qbuff, tri_surf, quad_surf);
        break;

      case Prism:
        insert_surf_pri(nod, ref_eidx[eidx], qbuff, tri_surf, quad_surf);
        break;

      case Hexa:
        insert_surf_hex(nod, ref_eidx[eidx], qbuff, quad_surf);
        break;

      default: break;
    }

    for(auto it = tri_surf.begin(); it != tri_surf.end(); ++it) {
      auto ft = search_tri.find(it->first);
      if(ft != search_tri.end() && !found_tri.count(ft->first)) {
        // now we also look for the same orientation
        T iv1 = mesh.pl.localize(it->second.v1);
        T iv2 = mesh.pl.localize(it->second.v2);
        T iv3 = mesh.pl.localize(it->second.v3);
        SF::Point ip1 = {mesh.xyz[iv1*3+0], mesh.xyz[iv1*3+1], mesh.xyz[iv1*3+2]};
        SF::Point ip2 = {mesh.xyz[iv2*3+0], mesh.xyz[iv2*3+1], mesh.xyz[iv2*3+2]};
        SF::Point ip3 = {mesh.xyz[iv3*3+0], mesh.xyz[iv3*3+1], mesh.xyz[iv3*3+2]};
        SF::Point e12 = normalize(ip2 - ip1);
        SF::Point e13 = normalize(ip3 - ip1);
        SF::Point in   = normalize(cross(e12, e13));

        T fv1 = mesh.pl.localize(ft->second.v1);
        T fv2 = mesh.pl.localize(ft->second.v2);
        T fv3 = mesh.pl.localize(ft->second.v3);
        SF::Point fp1 = {mesh.xyz[fv1*3+0], mesh.xyz[fv1*3+1], mesh.xyz[fv1*3+2]};
        SF::Point fp2 = {mesh.xyz[fv2*3+0], mesh.xyz[fv2*3+1], mesh.xyz[fv2*3+2]};
        SF::Point fp3 = {mesh.xyz[fv3*3+0], mesh.xyz[fv3*3+1], mesh.xyz[fv3*3+2]};
        e12 = normalize(fp2 - fp1);
        e13 = normalize(fp3 - fp1);
        SF::Point fn  = normalize(cross(e12, e13));

        if (inner_prod(in, fn) > 0.9) {
          tri_sele<T> found = ft->second; found.eidx = eidx;
          found_tri[ft->first] = found;
        }
      }
    }

    for(auto it = quad_surf.begin(); it != quad_surf.end(); ++it) {
      auto ft = search_quad.find(it->first);
      if(ft != search_quad.end() && !found_quad.count(ft->first)) {
        // now we also look for the same orientation
        T iv1 = mesh.pl.localize(it->second.v1);
        T iv2 = mesh.pl.localize(it->second.v2);
        T iv3 = mesh.pl.localize(it->second.v3);
        SF::Point ip1 = {mesh.xyz[iv1*3+0], mesh.xyz[iv1*3+1], mesh.xyz[iv1*3+2]};
        SF::Point ip2 = {mesh.xyz[iv2*3+0], mesh.xyz[iv2*3+1], mesh.xyz[iv2*3+2]};
        SF::Point ip3 = {mesh.xyz[iv3*3+0], mesh.xyz[iv3*3+1], mesh.xyz[iv3*3+2]};
        SF::Point e12 = normalize(ip2 - ip1);
        SF::Point e13 = normalize(ip3 - ip1);
        SF::Point in   = normalize(cross(e12, e13));

        T fv1 = mesh.pl.localize(ft->second.v1);
        T fv2 = mesh.pl.localize(ft->second.v2);
        T fv3 = mesh.pl.localize(ft->second.v3);
        SF::Point fp1 = {mesh.xyz[fv1*3+0], mesh.xyz[fv1*3+1], mesh.xyz[fv1*3+2]};
        SF::Point fp2 = {mesh.xyz[fv2*3+0], mesh.xyz[fv2*3+1], mesh.xyz[fv2*3+2]};
        SF::Point fp3 = {mesh.xyz[fv3*3+0], mesh.xyz[fv3*3+1], mesh.xyz[fv3*3+2]};
        e12 = normalize(fp2 - fp1);
        e13 = normalize(fp3 - fp1);
        SF::Point fn  = normalize(cross(e12, e13));

        if (inner_prod(in, fn) > 0.9) {
          quad_sele<T> found = ft->second; found.eidx = eidx;
          found_quad[ft->first] = found;
        }
      }
    }

    nod += mesh.dsp[eidx+1] - mesh.dsp[eidx];
  }
}

template<class T, class S> inline
void compute_surface(const meshdata<T,S> & mesh, const SF_nbr numbering,
                     const hashmap::unordered_set<T> & tags,
                     hashmap::unordered_map<triple<T>, tri_sele<T> >     & tri_surf,
                     hashmap::unordered_map<quadruple<T>, quad_sele<T> > & quad_surf)
{
  const T* con = mesh.con.data();
  const T* nbr = mesh.get_numbering(numbering).data();

  bool have_tags = tags.size() > 0;

  vector<T> nodvec(mesh.con.size()); T* nod = nodvec.data();
  for(size_t i=0; i<nodvec.size(); i++)
   nod[i] = nbr[con[i]];

  vector<T> qbuff(4); // buffer if quad nodes need to be sorted

  for(size_t eidx = 0; eidx < mesh.l_numelem; eidx++) {
    if(have_tags == false || tags.count(mesh.tag[eidx])) {
      switch(mesh.type[eidx]) {
        case Tri: {
          triple<T> tri;
          tri_sele<T> trie;
          insert_surf_tri(nod[0], nod[1], nod[2], eidx, tri, trie, tri_surf);
          break;
        }

        case Quad: {
          quadruple<T> qd;
          quad_sele<T> qde;
          insert_surf_quad(nod[0], nod[1], nod[2], nod[3],
                           eidx, qbuff, qd, qde, quad_surf);
          break;
        }

        case Tetra:
          insert_surf_tet(nod, eidx, tri_surf);
          break;

        case Pyramid:
          insert_surf_pyr(nod, eidx, qbuff, tri_surf, quad_surf);
          break;

        case Prism:
          insert_surf_pri(nod, eidx, qbuff, tri_surf, quad_surf);
          break;

        case Hexa:
          insert_surf_hex(nod, eidx, qbuff, quad_surf);
          break;

        default: break;
      }
    }

    nod += mesh.dsp[eidx+1] - mesh.dsp[eidx];
  }
}

template<class T, class S> inline
void compute_surface(const meshdata<T,S> & mesh, const SF_nbr numbering,
                     hashmap::unordered_map<triple<T>, tri_sele<T> >     & tri_surf,
                     hashmap::unordered_map<quadruple<T>, quad_sele<T> > & quad_surf)
{
  hashmap::unordered_set<T> tags;
  compute_surface(mesh, numbering, tags, tri_surf, quad_surf);
}

template<class V> inline
void get_hashmap_duplicates(const vector<V> & data, const MPI_Comm comm,
                            vector<bool> & is_dup)
{
  int size, rank;
  MPI_Comm_size(comm, &size); MPI_Comm_rank(comm, &rank);

  size_t dsize = data.size();
  is_dup.assign(dsize, false);

  vector<int> perm, dest(dsize);

  // set up destinations
  for(size_t i=0; i<dsize; i++)
    dest[i] = hashmap::hash_ops<V>::hash(data[i]) % size;

  commgraph<size_t> grph;
  grph.configure(dest, comm);
  size_t nrecv = sum(grph.rcnt);

  interval(perm, 0, dsize);
  binary_sort_copy(dest, perm);

  // fill send buffer and communicate
  vector<V> sbuff(dsize), rbuff(nrecv);
  for(size_t i=0; i<dsize; i++)
    sbuff[i] = data[perm[i]];
  MPI_Exchange(grph, sbuff, rbuff, comm);

  // now check duplicates. values occurring multiple times should go into dup, those occurring
  // only once go into not_dup
  hashmap::unordered_set<V> not_dup, dup;
  vector<bool> cdup(rbuff.size());

  for(const V & val : rbuff) {
    if(not_dup.count(val) == 0 && dup.count(val) == 0)
      not_dup.insert(val);
    else {
      not_dup.erase(val);
      dup.insert(val);
    }
  }

  for(size_t i=0; i<nrecv; i++) {
    V val = rbuff[i];
    bool d = dup.count(val);

    // just a debug check
    assert(d == (not_dup.count(val) == 0));

    cdup[i] = d;
  }

  // now we communicate cdup back
  grph.transpose();
  MPI_Exchange(grph, cdup, is_dup, comm);

  // permute is_dup back to input permutation
  binary_sort_copy(perm, is_dup);
}

template<class K, class V> inline
void remove_duplicates(hashmap::unordered_map<K,V> & map, const MPI_Comm comm)
{
  vector<K>    check_vec (map.size());
  vector<bool> is_dup    (map.size());

  size_t idx=0;
  for(const auto & v : map) check_vec[idx++] = v.first;

  get_hashmap_duplicates(check_vec, comm, is_dup);

  for(size_t i=0; i<is_dup.size(); i++)
    if(is_dup[i])
      map.erase(check_vec[i]);
}
template<class K> inline
void remove_duplicates(hashmap::unordered_set<K> & set, const MPI_Comm comm)
{
  vector<K>    check_vec (set.size());
  vector<bool> is_dup    (set.size());

  check_vec.assign(set.begin(), set.end());
  get_hashmap_duplicates(check_vec, comm, is_dup);

  for(size_t i=0; i<is_dup.size(); i++)
    if(is_dup[i])
      set.erase(check_vec[i]);
}

template<class T> inline
void remove_parallel_duplicates(hashmap::unordered_map<triple<T>,tri_sele<T>>      & tri_surf,
                                hashmap::unordered_map<quadruple<T>,quad_sele<T> > & quad_surf,
                                MPI_Comm comm)
{
  long int g_num_tri  = tri_surf.size();
  long int g_num_quad = quad_surf.size();
  MPI_Allreduce(MPI_IN_PLACE, &g_num_tri,  1, MPI_LONG, MPI_SUM, comm);
  MPI_Allreduce(MPI_IN_PLACE, &g_num_quad, 1, MPI_LONG, MPI_SUM, comm);

  if(g_num_tri)
    remove_duplicates(tri_surf, comm);
  if(g_num_quad)
    remove_duplicates(quad_surf, comm);
}


template<class T, class S> inline
void convert_surface_mesh(hashmap::unordered_map<triple<T>,tri_sele<T>>      & tri_surf,
                          hashmap::unordered_map<quadruple<T>,quad_sele<T> > & quad_surf,
                          meshdata<T,S> & surfmesh,
                          vector<T> & elem_orig)
{
  MPI_Comm comm = surfmesh.comm;

  long int g_num_tri  = tri_surf.size();
  long int g_num_quad = quad_surf.size();
  MPI_Allreduce(MPI_IN_PLACE, &g_num_tri,  1, MPI_LONG, MPI_SUM, comm);
  MPI_Allreduce(MPI_IN_PLACE, &g_num_quad, 1, MPI_LONG, MPI_SUM, comm);

  surfmesh.g_numelem = g_num_tri + g_num_quad;
  surfmesh.l_numelem = tri_surf.size() + quad_surf.size();

  vector<T> cnt(surfmesh.l_numelem);
  surfmesh.type.resize(surfmesh.l_numelem);
  surfmesh.con.resize(tri_surf.size() * 3 + quad_surf.size() * 4);
  elem_orig.resize(surfmesh.l_numelem);

  size_t idx = 0, cidx = 0;
  for(const auto & v : tri_surf) {
    cnt[idx] = 3;
    surfmesh.type[idx] = Tri;
    surfmesh.con[cidx + 0] = v.second.v1;
    surfmesh.con[cidx + 1] = v.second.v2;
    surfmesh.con[cidx + 2] = v.second.v3;

    elem_orig[idx] = v.second.eidx;

    idx  += 1;
    cidx += 3;
  }

  for(const auto & v : quad_surf) {
    cnt[idx] = 4;
    surfmesh.type[idx] = Quad;
    surfmesh.con[cidx + 0] = v.second.v1;
    surfmesh.con[cidx + 1] = v.second.v2;
    surfmesh.con[cidx + 2] = v.second.v3;
    surfmesh.con[cidx + 3] = v.second.v4;

    idx  += 1;
    cidx += 4;
  }

  dsp_from_cnt(cnt, surfmesh.dsp);
}

template<class T, class S> inline
void convert_mesh_surface(const meshdata<T,S> & surfmesh,
                          hashmap::unordered_map<triple<T>,tri_sele<T>>     & tri_surf,
                          hashmap::unordered_map<quadruple<T>,quad_sele<T>> & quad_surf)
{
  vector<T> qbuff(4); // buffer if quad nodes need to be sorted
  quadruple<T> qd;  quad_sele<T> qde;
  triple<T>    tri; tri_sele<T>  trie;

  for(size_t eidx=0; eidx<surfmesh.l_numelem; eidx++) {
    int esize = surfmesh.dsp[eidx+1] - surfmesh.dsp[eidx];

    const T* nod = surfmesh.con.data() + surfmesh.dsp[eidx];
    if(esize == 4) {
      insert_surf_quad(nod[0], nod[1], nod[2], nod[3], eidx, qbuff, qd, qde, quad_surf);
    }
    else {
      insert_surf_tri(nod[0], nod[1], nod[2], eidx, tri, trie, tri_surf);
    }
  }
}

template<class T, class S> inline
void compute_surface_mesh(const meshdata<T,S> & mesh, const SF_nbr numbering,
                          const hashmap::unordered_set<T> & tags,
                          meshdata<T,S> & surfmesh)
{
  hashmap::unordered_map<triple<T>, tri_sele<T> >      tri_surf;
  hashmap::unordered_map<quadruple<T>, quad_sele<T> >  quad_surf;
  vector<T> elem_orig; // local element indices of "mesh"

  compute_surface(mesh, numbering, tags, tri_surf, quad_surf);

  surfmesh.comm = mesh.comm;
  remove_parallel_duplicates(tri_surf, quad_surf, surfmesh.comm);
  convert_surface_mesh(tri_surf, quad_surf, surfmesh, elem_orig);

  surfmesh.tag.resize(size_t(surfmesh.l_numelem));
  surfmesh.fib.resize(size_t(surfmesh.l_numelem * 3));

  if(mesh.she.size())
    surfmesh.she.resize(surfmesh.l_numelem * 3);

  for(size_t eidx = 0; eidx < elem_orig.size(); eidx++) {
    T orig = elem_orig[eidx];

    surfmesh.tag[eidx]     = mesh.tag[orig];
    surfmesh.fib[eidx*3+0] = mesh.fib[orig*3+0];
    surfmesh.fib[eidx*3+1] = mesh.fib[orig*3+1];
    surfmesh.fib[eidx*3+2] = mesh.fib[orig*3+2];

    if(mesh.she.size()) {
      surfmesh.she[eidx*3+0] = mesh.she[orig*3+0];
      surfmesh.she[eidx*3+1] = mesh.she[orig*3+1];
      surfmesh.she[eidx*3+2] = mesh.she[orig*3+2];
    }
  }
}

template<class T, class S> inline
void compute_surface_mesh(const meshdata<T,S> & mesh, const SF_nbr numbering,
                          meshdata<T,S> & surfmesh)
{
  hashmap::unordered_set<T> tags;
  compute_surface_mesh(mesh, numbering, tags, surfmesh);
}

template<class T, class S> inline
void read_surface_mesh(const meshdata<T,S> & mesh,
                       meshdata<T,S> & surface,
                       std::string filename)
{
  int rank; MPI_Comm_rank(mesh.comm, &rank);
  FILE *fd = NULL;

  size_t numele   = 0;
  size_t numcon   = 0;
  bool   twoFib   = false;
  bool   read_bin = false;

  if(rank == 0)
  {
    fd = fopen(filename.c_str(), "r");
    if(fd)
      SF::read_headers(fd, NULL, read_bin, numele, twoFib);
  }
  MPI_Bcast(&numele, sizeof(size_t), MPI_BYTE, 0, mesh.comm);
  surface.g_numelem = numele;

  // We read the whole surface at once, this could be rewritten into multiple reads to
  // save memory
  if(rank == 0) {
    SF::read_elem_block(fd, read_bin, 0, numele, surface);

    numcon = surface.con.size();

    if(numele != surface.l_numelem) {
      std::cerr << "Error: Incomplete surface file! Aborting!" << std::endl;
      fclose(fd);
      exit(1);
    }
  }

  surface.l_numelem = numele;

  surface.tag.resize(surface.l_numelem);
  MPI_Bcast(surface.tag.data(), surface.tag.size()*sizeof(T), MPI_BYTE, 0, mesh.comm);

  surface.dsp.resize(surface.l_numelem + 1);
  MPI_Bcast(surface.dsp.data(), surface.dsp.size()*sizeof(T), MPI_BYTE, 0, mesh.comm);

  surface.type.resize(surface.l_numelem);
  MPI_Bcast(surface.type.data(), surface.type.size()*sizeof(elem_t), MPI_BYTE, 0, mesh.comm);

  vector<T> & ref_eidx = surface.register_numbering(NBR_ELEM_REF);
  ref_eidx.resize(surface.l_numelem);

  MPI_Bcast(ref_eidx.data(), surface.l_numelem*sizeof(T), MPI_BYTE, 0, mesh.comm);

  MPI_Bcast(&numcon, sizeof(size_t), MPI_BYTE, 0, mesh.comm);
  surface.con.resize(numcon);
  MPI_Bcast(surface.con.data(), surface.con.size()*sizeof(T), MPI_BYTE, 0, mesh.comm);

  // now we have to restrict the surface elements to those in the local partition
  {
    size_t initial_gnumelem = surface.g_numelem;

    hashmap::unordered_map<triple<T>, tri_sele<T>>     search_tri,  found_tri;
    hashmap::unordered_map<quadruple<T>, quad_sele<T>> search_quad, found_quad;

    convert_mesh_surface(surface, search_tri, search_quad);
    search_for_surface  (mesh, SF::NBR_REF, search_tri, search_quad, found_tri, found_quad);

    vector<T> elem_orig;
    convert_surface_mesh(found_tri, found_quad, surface, elem_orig);
    mesh.pl.localize(surface.con);

    surface.tag.resize(size_t(surface.l_numelem));
    surface.fib.resize(size_t(surface.l_numelem * 3));

    if(mesh.she.size())
      surface.she.resize(surface.l_numelem * 3);

    for(size_t eidx = 0; eidx < elem_orig.size(); eidx++) {
      T orig = elem_orig[eidx];

      surface.tag[eidx]     = mesh.tag[orig];
      surface.fib[eidx*3+0] = mesh.fib[orig*3+0];
      surface.fib[eidx*3+1] = mesh.fib[orig*3+1];
      surface.fib[eidx*3+2] = mesh.fib[orig*3+2];

      if(mesh.she.size()) {
        surface.she[eidx*3+0] = mesh.she[orig*3+0];
        surface.she[eidx*3+1] = mesh.she[orig*3+1];
        surface.she[eidx*3+2] = mesh.she[orig*3+2];
      }
    }

    // we want to remember to original gnumelems, so that the next consistency check is
    // actually useful.
    surface.g_numelem = initial_gnumelem;
  }

  {
    long int numele_check = surface.l_numelem;
    MPI_Allreduce(MPI_IN_PLACE, &numele_check, 1, MPI_LONG, MPI_SUM, mesh.comm);
    if(numele_check != surface.g_numelem) {
      if(rank == 0)
        fprintf(stderr, "ERROR: Bad partitioning of surface %s!"
                        " Global elem sum should be %ld, but is %ld!\n\n",
                filename.c_str(), surface.g_numelem, numele_check);
    }
  }

  // copy indexing from parent mesh
  SF::vector<T> & nod_ref   = surface.register_numbering(SF::NBR_REF);
  SF::vector<T> & nod_petsc = surface.register_numbering(SF::NBR_PETSC);
  nod_ref   = mesh.get_numbering(SF::NBR_REF);
  nod_petsc = mesh.get_numbering(SF::NBR_PETSC);

  // copy coords from parent mesh
  surface.l_numpts = mesh.l_numpts;
  surface.g_numpts = mesh.g_numpts;
  surface.xyz = mesh.xyz;
}

template<class T>
inline void restrict_to_set(vector<T> & v, const hashmap::unordered_set<T> & set)
{
  size_t widx = 0;

  for(size_t i=0; i<v.size(); i++)
    if(set.count(v[i]))
      v[widx++] = v[i];

  v.resize(widx);
}
template<class T>
inline void restrict_to_set(vector<T> & v, const vector<T> & vset)
{
  hashmap::unordered_set<T> set;
  set.insert(vset.begin(), vset.end());

  restrict_to_set(v, set);
}

}

#endif