dist_grid.c

Go to the documentation of this file.

 
/*
 * Keywords:
 * Maintainer: Moritz Hanke <hanke@dkrz.de>
 * URL: https://dkrz-sw.gitlab-pages.dkrz.de/yac/
 *
 * This file is part of YAC.
 *
 * Redistribution and use in source and binary forms, with or without
 * modification, are  permitted provided that the following conditions are
 * met:
 *
 * Redistributions of source code must retain the above copyright notice,
 * this list of conditions and the following disclaimer.
 *
 * Redistributions in binary form must reproduce the above copyright
 * notice, this list of conditions and the following disclaimer in the
 * documentation and/or other materials provided with the distribution.
 *
 * Neither the name of the DKRZ GmbH nor the names of its contributors
 * may be used to endorse or promote products derived from this software
 * without specific prior written permission.
 *
 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS "AS
 * IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
 * TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A
 * PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT OWNER
 * OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL,
 * EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO,
 * PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR
 * PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF
 * LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING
 * NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS
 * SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
 */
 
#ifdef HAVE_CONFIG_H
// Get the definition of the 'restrict' keyword.
#include "config.h"
#endif
 
#include <stdlib.h>
#include <string.h>
#include <assert.h>
#include <float.h>
#include <yaxt.h>
 
#include "grid.h"
#include "dist_grid.h"
#include "geometry.h"
#include "yac_mpi.h"
#include "utils.h"
#include "sphere_part.h"
#include "ensure_array_size.h"
#include "interp_grid.h"
 
struct single_remote_point {
  yac_int global_id;
  struct remote_point_info data;
};
 
struct single_remote_point_reorder {
  struct single_remote_point data;
  size_t reorder_idx;
};
 
struct remote_point_reorder {
  struct remote_point data;
  size_t reorder_idx;
};
 
// reorder_idx has to be first
struct global_vertex_reorder {
  size_t reorder_idx;
  yac_int global_id;
  double coord[3];
  struct remote_point_infos owners;
};
 
// reorder_idx has to be first
struct global_edge_reorder {
  size_t reorder_idx;
  yac_int global_id;
  enum yac_edge_type edge_type;
  yac_int edge_to_vertex[2];
  struct remote_point_infos owners;
};
 
struct id_pos_coord {
  yac_int global_id;
  uint64_t orig_pos;
  double coord[3];
};
 
struct id_pos {
  yac_int global_id;
  uint64_t orig_pos;
};
 
struct n_ids_reorder {
 
  int count;
  yac_int * ids;
  size_t reorder_idx;
  yac_int global_id;
};
 
struct id_reorder_coord {
  yac_int global_id;
  size_t reorder_idx;
  double coord[3];
};
 
struct temp_field_data {
  int ** masks;
  size_t * masks_array_sizes;
  size_t masks_count;
  coordinate_pointer * coordinates;
  size_t * coordinates_array_sizes;
  size_t coordinates_count;
};
 
// warning: when changing this, ensure that struct const_basic_grid_data is
// changed accordingly
struct dist_grid {
  coordinate_pointer vertex_coordinates;
  yac_int * cell_ids;
  yac_int * vertex_ids;
  yac_int * edge_ids;
  int * num_vertices_per_cell;
  size_t * cell_to_vertex;
  size_t * cell_to_vertex_offsets;
  size_t * cell_to_edge;
  size_t * cell_to_edge_offsets;
  size_t_2_pointer edge_to_vertex;
  struct bounding_circle * cell_bnd_circles;
  enum yac_edge_type * edge_type;
  struct remote_point_infos * cell_owners;   // original cell owner(s)
  struct remote_point_infos * vertex_owners; // original vertex owner(s)
  struct remote_point_infos * edge_owners;   // original edge owner(s)
  size_t num_total_cells;
  size_t num_total_vertices;
  size_t num_total_edges;
  size_t num_cells;
  size_t num_vertices;
  size_t num_edges;
  int * core_cell_mask;
  int * core_vertex_mask;
  int * core_edge_mask;
  yac_int * sorted_cell_ids;
  yac_int * sorted_vertex_ids;
  yac_int * sorted_edge_ids;
  size_t * sorted_cell_reorder_idx;
  size_t * sorted_vertex_reorder_idx;
  size_t * sorted_edge_reorder_idx;
  struct yac_field_data_set field_data;
  MPI_Comm comm;
};
 
struct dist_grid_pair {
  struct dist_grid dist_grid[2];
  int grid_ids[2];
  struct proc_sphere_part_node * proc_sphere_part;
  struct point_sphere_part_search * vertex_sphere_part[2];
  struct bnd_sphere_part_search * cell_sphere_part[2];
  MPI_Comm comm;
};
 
struct nnn_search_result {
  struct single_remote_point point_info;
  double cos_angle;
};
 
struct missing_edge_neighbour {
  struct {
    size_t local_id;
    yac_int global_id;
  } cell, edge;
  size_t neigh_idx;
};
 
// looks up positions of ids in an array of sorted ids
static void id2idx(
  yac_int * ids, size_t * idx, size_t num_ids,
  yac_int * ref_sorted_ids, size_t num_sorted_ids) {
 
  size_t * reorder = xmalloc(num_ids * sizeof(*reorder));
  for (size_t i = 0; i < num_ids; ++i) reorder[i] = i;
 
  yac_quicksort_index_yac_int_size_t(ids, num_ids, reorder);
 
  for (size_t i = 0, j = 0; i < num_ids; ++i) {
 
    yac_int curr_id = ids[i];
    while ((j < num_sorted_ids) && (ref_sorted_ids[j] < curr_id)) ++j;
    YAC_ASSERT(
      (j < num_sorted_ids) && (ref_sorted_ids[j] == curr_id),
      "ERROR(id2idx): id not found")
    idx[reorder[i]] = j;
  }
 
  free(reorder);
}
 
static int * determine_core_vertex_mask(
  struct proc_sphere_part_node * proc_sphere_part, int comm_rank,
  coordinate_pointer vertex_coordinates, size_t count) {
 
  int * core_mask = xcalloc(count, sizeof(*core_mask));
 
  yac_proc_sphere_part_do_point_search(
    proc_sphere_part, vertex_coordinates, count, core_mask);
 
  for (size_t i = 0; i < count; ++i)
    core_mask[i] = core_mask[i] == comm_rank;
 
  return core_mask;
}
 
static size_t get_cell_reference_vertex(
  struct dist_grid * grid, size_t cell_idx) {
 
  int num_vertices = grid->num_vertices_per_cell[cell_idx];
  if (num_vertices == 0) return SIZE_MAX;
  yac_int * grid_vertex_ids = grid->vertex_ids;
  size_t * vertices =
    grid->cell_to_vertex + grid->cell_to_vertex_offsets[cell_idx];
  // get the cell corner with the smallest global id
  size_t min_idx = vertices[0];
  yac_int min_global_id = grid_vertex_ids[min_idx];
  for (int j = 1; j < num_vertices; ++j) {
    size_t curr_vertex = vertices[j];
    yac_int curr_global_id = grid_vertex_ids[curr_vertex];
    if (min_global_id > curr_global_id) {
      min_global_id = curr_global_id;
      min_idx = curr_vertex;
    }
  }
  return min_idx;
}
 
// determines the core mask for the cells in the distributed grid
// the core mask value is equal to that of the core mask value of the
// cell vertex with the lowest global id
static int * determine_cell_core_mask(
  int const * vertex_core_mask, struct dist_grid * dist_grid, int is_root) {
 
  size_t num_cells = dist_grid->num_cells;
  int * core_cell_mask = xmalloc(num_cells * sizeof(core_cell_mask));
 
  //-------------------------
  // determine cell core mask
  //-------------------------
  for (size_t i = 0; i < num_cells; ++i) {
    size_t ref_vertex = get_cell_reference_vertex(dist_grid, i);
    core_cell_mask[i] =
      (ref_vertex != SIZE_MAX)?(vertex_core_mask[ref_vertex]):is_root;
  }
  return core_cell_mask;
}
 
static inline size_t get_edge_reference_vertex(
  struct dist_grid * grid, size_t edge_idx) {
 
  yac_int * vertex_ids = grid->vertex_ids;
  size_t * edge_vertices = &(grid->edge_to_vertex[edge_idx][0]);
  // get the edge corner with the smallest global id
  return edge_vertices[
          (vertex_ids[edge_vertices[0]] > vertex_ids[edge_vertices[1]])?1:0];
}
 
// determines the core mask for the edges in the distributed grid
// the core mask value is equal to that of the core mask value of the
// edge vertex with the lowest global id
static int * determine_edge_core_mask(
  int const * vertex_core_mask, struct dist_grid * dist_grid) {
 
  size_t num_edges = dist_grid->num_edges;
  int * core_edge_mask = xmalloc(num_edges * sizeof(core_edge_mask));
 
  //-------------------------
  // determine edge core mask
  //-------------------------
  for (size_t i = 0; i < num_edges; ++i)
    core_edge_mask[i] =
      vertex_core_mask[get_edge_reference_vertex(dist_grid, i)];
  return core_edge_mask;
}
 
static MPI_Datatype yac_get_single_remote_point_mpi_datatype(MPI_Comm comm) {
 
  struct single_remote_point dummy;
  MPI_Datatype single_id_owner_dt;
  int array_of_blocklengths[] = {1, 1, 1};
  const MPI_Aint array_of_displacements[] =
    {(MPI_Aint)((void*)&(dummy.global_id) - (void*)&dummy),
     (MPI_Aint)((void*)&(dummy.data.rank) - (void*)&dummy),
     (MPI_Aint)((void*)&(dummy.data.orig_pos) - (void*)&dummy)};
  const MPI_Datatype array_of_types[] =
    {Xt_int_dt, MPI_INT, MPI_UINT64_T};
  yac_mpi_call(
    MPI_Type_create_struct(3, array_of_blocklengths, array_of_displacements,
                           array_of_types, &single_id_owner_dt), comm);
  return yac_create_resized(single_id_owner_dt, sizeof(dummy), comm);
}
 
static MPI_Datatype yac_get_n_single_remote_point_reorder_mpi_datatype(
  int n, MPI_Comm comm) {
 
  MPI_Datatype dt_single_remote_point_ =
    yac_get_single_remote_point_mpi_datatype(comm);
  MPI_Datatype dt_single_remote_point__ =
    yac_create_resized(
      dt_single_remote_point_, sizeof(struct single_remote_point_reorder), comm);
  MPI_Datatype dt_single_remote_point;
  yac_mpi_call(
    MPI_Type_contiguous(
      n, dt_single_remote_point__, &dt_single_remote_point), comm);
  yac_mpi_call(MPI_Type_free(&dt_single_remote_point__), comm);
  yac_mpi_call(MPI_Type_commit(&dt_single_remote_point), comm);
  return dt_single_remote_point;
}
 
static int compare_single_remote_point_global_id (
  const void * a, const void * b) {
 
  return (((const struct single_remote_point *)a)->global_id >
          ((const struct single_remote_point *)b)->global_id) -
         (((const struct single_remote_point *)a)->global_id <
          ((const struct single_remote_point *)b)->global_id);
}
 
static size_t determine_local_count(size_t count, int * owner) {
 
  size_t local_count = 0;
  for (size_t i = 0; i < count; ++i) if (owner[i]) ++local_count;
  return local_count;
}
 
static MPI_Datatype yac_get_id_pos_mpi_datatype(MPI_Comm comm) {
 
  struct id_pos dummy;
  MPI_Datatype id_pos_dt;
  int array_of_blocklengths[] = {1,1};
  const MPI_Aint array_of_displacements[] =
    {(MPI_Aint)((void*)&(dummy.global_id) - (void*)&dummy),
     (MPI_Aint)((void*)&(dummy.orig_pos) - (void*)&dummy)};
  const MPI_Datatype array_of_types[] = {Xt_int_dt, MPI_UINT64_T};
  yac_mpi_call(
    MPI_Type_create_struct(
      2, array_of_blocklengths, array_of_displacements,
      array_of_types, &id_pos_dt), comm);
  return yac_create_resized(id_pos_dt, sizeof(dummy), comm);
}
 
MPI_Datatype yac_get_remote_point_info_mpi_datatype(MPI_Comm comm) {
 
  struct remote_point_info dummy;
  MPI_Datatype remote_point_info_dt;
  int array_of_blocklengths[] = {1, 1};
  const MPI_Aint array_of_displacements[] =
    {(MPI_Aint)((void*)&(dummy.rank) - (void*)&dummy),
     (MPI_Aint)((void*)&(dummy.orig_pos) - (void*)&dummy)};
  const MPI_Datatype array_of_types[] =
    {MPI_INT, MPI_UINT64_T};
  yac_mpi_call(
    MPI_Type_create_struct(2, array_of_blocklengths, array_of_displacements,
                           array_of_types, &remote_point_info_dt), comm);
  return yac_create_resized(remote_point_info_dt, sizeof(dummy), comm);
}
 
static struct remote_point_infos copy_remote_point_infos(
  struct remote_point_infos point_infos) {
 
  int count = point_infos.count;
  if (count > 1) {
    struct remote_point_info * temp = xmalloc((size_t)count * sizeof(*temp));
    memcpy(temp, point_infos.data.multi, (size_t)count * sizeof(*temp));
    point_infos.data.multi = temp;
  }
  return point_infos;
}
 
static struct remote_point_infos * generate_remote_point_infos(
  yac_int * sorted_ids, size_t * reorder_idx, size_t count,
  struct remote_points * points) {
 
  struct remote_point_infos * point_infos = xmalloc(count * sizeof(*point_infos));
 
  size_t id_owner_count = points->count;
  struct remote_point * points_ = points->data;
 
  for (size_t i = 0, j = 0; i < count; ++i) {
 
    yac_int curr_id = sorted_ids[i];
 
    while ((j < id_owner_count) && (points_[j].global_id < curr_id)) ++j;
 
    if ((j >= id_owner_count) || (points_[j].global_id != curr_id))
      point_infos[reorder_idx[i]].count = -1;
    else
      point_infos[reorder_idx[i]] = copy_remote_point_infos(points_[j].data);
  }
 
  return point_infos;
}
 
static void lookup_remote_points(
  yac_int * requested_ids, struct remote_point_infos ** results,
  size_t * reorder_buffer, size_t request_count,
  struct remote_points * local_remote_data) {
 
  struct remote_point * local_remote_points = local_remote_data->data;
  size_t local_count = local_remote_data->count;
  static struct remote_point_infos dummy_data =
    {.count = 0, .data.single = {.rank = INT_MAX, .orig_pos = UINT64_MAX}};
 
  for (size_t j = 0; j < request_count; ++j) reorder_buffer[j] = j;
 
  // sort ids for easier owner lookup
  yac_quicksort_index_yac_int_size_t(
    requested_ids, request_count, reorder_buffer);
 
  for (size_t i = 0, j = 0; i < request_count; ++i) {
    yac_int curr_id = requested_ids[i];
    while ((j < local_count) && (curr_id > local_remote_points[j].global_id))
      ++j;
    results[reorder_buffer[i]] =
      ((j < local_count) && (curr_id == local_remote_points[j].global_id))?
        &(local_remote_points[j].data):&dummy_data;
  }
}
 
static int yac_remote_point_infos_get_pack_size(
  struct remote_point_infos const * infos, MPI_Datatype point_info_dt,
  MPI_Comm comm) {
 
  int pack_size_count,
      pack_size_data;
 
  yac_mpi_call(MPI_Pack_size(1, MPI_INT, comm, &pack_size_count), comm);
  yac_mpi_call(
    MPI_Pack_size(infos->count, point_info_dt, comm, &pack_size_data), comm);
 
  return pack_size_count + pack_size_data;
}
 
static void yac_remote_point_infos_get_pack_sizes(
  struct remote_point_infos ** infos, size_t count, int * pack_sizes,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  int pack_size_count;
 
  yac_mpi_call(MPI_Pack_size(1, MPI_INT, comm, &pack_size_count), comm);
 
  for (size_t i = 0; i < count; ++i) {
    yac_mpi_call(
      MPI_Pack_size(
        infos[i]->count, point_info_dt, comm, pack_sizes + i), comm);
    pack_sizes[i] += pack_size_count;
  }
}
 
static void yac_remote_point_infos_pack(
  struct remote_point_infos ** infos, int * pack_sizes, size_t count,
  void * buffer, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  for (size_t i = 0; i < count; ++i) {
 
    int curr_count = infos[i]->count;
    struct remote_point_info * curr_point_infos =
      (curr_count == 1)?
        (&(infos[i]->data.single)):(infos[i]->data.multi);
 
    int position = 0;
    int buffer_size = pack_sizes[i];
    yac_mpi_call(
        MPI_Pack(&curr_count, 1, MPI_INT, buffer,
                 buffer_size, &position, comm), comm);
    yac_mpi_call(
        MPI_Pack(curr_point_infos, curr_count, point_info_dt, buffer,
                 buffer_size, &position, comm), comm);
    buffer = (void*)(((unsigned char*)buffer) + buffer_size);
  }
}
 
static void yac_remote_point_infos_unpack(
  void * buffer, int buffer_size, int * position,
  struct remote_point_infos * infos, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  int count;
  yac_mpi_call(
    MPI_Unpack(buffer, buffer_size, position, &count, 1, MPI_INT, comm), comm);
 
  infos->count = count;
 
  if (count > 0) {
    struct remote_point_info * point_infos;
    if (count == 1)
      point_infos = &(infos->data.single);
    else
      point_infos =
        (infos->data.multi = xmalloc((size_t)count * sizeof(*point_infos)));
 
    yac_mpi_call(
      MPI_Unpack(buffer, buffer_size, position, point_infos, count,
                 point_info_dt, comm), comm);
  } else {
    infos->data.single.rank = INT_MAX;
    infos->data.single.orig_pos = UINT64_MAX;
  }
}
 
static void yac_remote_point_infos_unpack_info_buffer(
  void * buffer, int buffer_size, int * position,
  struct remote_point_info * info_buffer, size_t * info_buffer_position,
  struct remote_point_infos * infos, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  int count;
  yac_mpi_call(
    MPI_Unpack(buffer, buffer_size, position, &count, 1, MPI_INT, comm), comm);
 
  infos->count = count;
 
  struct remote_point_info * point_infos;
  if (count == 1)
    point_infos = &(infos->data.single);
  else {
    point_infos =
      (infos->data.multi = info_buffer + *info_buffer_position);
    *info_buffer_position += count;
  }
 
  yac_mpi_call(
    MPI_Unpack(buffer, buffer_size, position, point_infos, count,
               point_info_dt, comm), comm);
}
 
int yac_remote_point_get_pack_size(
  struct remote_point * point, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  int pack_size_id,
      pack_size_remote_point_infos;
 
  yac_mpi_call(MPI_Pack_size(1, Xt_int_dt, comm, &pack_size_id), comm);
  pack_size_remote_point_infos =
    yac_remote_point_infos_get_pack_size(&(point->data), point_info_dt, comm);
 
  return pack_size_id + pack_size_remote_point_infos;
}
 
static void yac_remote_point_infos_pack_single(
  struct remote_point_infos const * infos, void * buffer, int buffer_size,
  int * position, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  int count = infos->count;
 
  struct remote_point_info const * info =
    (count == 1)?(&(infos->data.single)):(infos->data.multi);
 
  yac_mpi_call(
      MPI_Pack(&count, 1, MPI_INT, buffer, buffer_size, position, comm), comm);
  yac_mpi_call(
      MPI_Pack(info, count, point_info_dt, buffer, buffer_size, position, comm),
      comm);
}
 
void yac_remote_point_pack(
  struct remote_point * point, void * buffer, int buffer_size, int * position,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  yac_mpi_call(
      MPI_Pack(&(point->global_id), 1, Xt_int_dt, buffer,
               buffer_size, position, comm), comm);
 
  yac_remote_point_infos_pack_single(
    &(point->data), buffer, buffer_size, position, point_info_dt, comm);
}
 
void yac_remote_point_unpack(
  void * buffer, int buffer_size, int * position, struct remote_point * point,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  yac_mpi_call(
    MPI_Unpack(
      buffer, buffer_size, position, &(point->global_id), 1, Xt_int_dt, comm),
    comm);
  yac_remote_point_infos_unpack(
    buffer, buffer_size, position, &(point->data), point_info_dt, comm);
}
 
void yac_remote_point_unpack_info_buffer(
  void * buffer, int buffer_size, int * position,
  struct remote_point_info * info_buffer, size_t * info_buffer_position,
  struct remote_point * point, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  yac_mpi_call(
    MPI_Unpack(
      buffer, buffer_size, position, &(point->global_id), 1, Xt_int_dt, comm),
    comm);
  yac_remote_point_infos_unpack_info_buffer(
    buffer, buffer_size, position, info_buffer, info_buffer_position,
    &(point->data), point_info_dt, comm);
}
 
int yac_remote_points_get_pack_size(
  struct remote_points * points, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  size_t count = points->count;
  struct remote_point * points_data = points->data;
 
  int count_pack_size,
      remote_points_pack_size;
 
  yac_mpi_call(MPI_Pack_size(2, MPI_UINT64_T, comm, &count_pack_size), comm);
  remote_points_pack_size = 0;
  for (size_t i = 0; i < count; ++i)
    remote_points_pack_size +=
      yac_remote_point_get_pack_size(points_data + i, point_info_dt, comm);
 
  return count_pack_size + remote_points_pack_size;
}
 
void yac_remote_points_pack(
  struct remote_points * points, void * buffer, int buffer_size, int * position,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  size_t count = points->count;
  struct remote_point * point_data = points->data;
  uint64_t counts[2] = {(uint64_t)count, 0};
  for (size_t i = 0; i < count; ++i)
    if (point_data[i].data.count > 1)
      counts[1] += (uint64_t)(point_data[i].data.count);
 
  yac_mpi_call(
      MPI_Pack(counts, 2, MPI_UINT64_T, buffer,
               buffer_size, position, comm), comm);
  for (size_t i = 0; i < count; ++i)
    yac_remote_point_pack(
      point_data + i, buffer, buffer_size, position, point_info_dt, comm);
}
 
void yac_remote_points_unpack(
  void * buffer, int buffer_size, int * position,
  struct remote_points ** points, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  uint64_t counts[2];
 
  yac_mpi_call(
    MPI_Unpack(
      buffer, buffer_size, position, counts, 2, MPI_UINT64_T, comm), comm);
 
  *points = xmalloc(((size_t)counts[1]) * sizeof(struct remote_point_infos) +
                    sizeof(**points));
 
  size_t count = ((*points)->count = (size_t)(counts[0]));
  struct remote_point * point_data =
    ((*points)->data = xmalloc(count * sizeof(*((*points)->data))));
  struct remote_point_info * remote_point_info_buffer =
    &((*points)->buffer[0]);
 
  for (size_t i = 0, offset = 0; i < count; ++i) {
 
    yac_remote_point_unpack_info_buffer(
      buffer, buffer_size, position, remote_point_info_buffer, &offset,
      point_data + i, point_info_dt, comm);
  }
}
 
static void generate_dist_grid_remote_point_infos(
  struct dist_grid * dist_grid, struct proc_sphere_part_node * proc_sphere_part,
  struct remote_points * dist_cell_owner,
  struct remote_points * dist_vertex_owner,
  struct remote_points * dist_edge_owner) {
 
  MPI_Comm comm = dist_grid->comm;
 
  int comm_size;
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  dist_grid->cell_owners =
    generate_remote_point_infos(
      dist_grid->sorted_cell_ids, dist_grid->sorted_cell_reorder_idx,
      dist_grid->num_cells, dist_cell_owner);
  dist_grid->vertex_owners =
    generate_remote_point_infos(
      dist_grid->sorted_vertex_ids, dist_grid->sorted_vertex_reorder_idx,
      dist_grid->num_vertices, dist_vertex_owner);
  dist_grid->edge_owners =
    generate_remote_point_infos(
      dist_grid->sorted_edge_ids, dist_grid->sorted_edge_reorder_idx,
      dist_grid->num_edges, dist_edge_owner);
 
  size_t num_missing_vertices = 0;
  size_t num_missing_edges = 0;
 
  for (size_t i = 0; i < dist_grid->num_cells; ++i)
    YAC_ASSERT(
      dist_grid->cell_owners[i].count != -1,
      "ERROR(generate_dist_grid_remote_point_infos):"
      "no owner information for a cell available")
  for (size_t i = 0; i < dist_grid->num_vertices; ++i)
    if (dist_grid->vertex_owners[i].count == -1) ++num_missing_vertices;
  for (size_t i = 0; i < dist_grid->num_edges; ++i)
    if (dist_grid->edge_owners[i].count == -1) ++num_missing_edges;
 
  size_t total_num_missing = num_missing_vertices + num_missing_edges;
  int * ranks_buffer = xmalloc(total_num_missing * sizeof(*ranks_buffer));
  int * vertex_ranks = ranks_buffer;
  int * edge_ranks = ranks_buffer + num_missing_vertices;
  size_t * size_t_buffer =
    xmalloc(3 * total_num_missing * sizeof(*size_t_buffer));
  size_t * vertex_reorder = size_t_buffer;
  size_t * edge_reorder = size_t_buffer + num_missing_vertices;
  size_t * coord_indices = size_t_buffer + total_num_missing;
  size_t * ranks_idxs = coord_indices;
  size_t * cve_reorder = size_t_buffer + 2 * total_num_missing;
 
  size_t offset = 0;
  for (size_t i = 0, j = 0; i < dist_grid->num_vertices; ++i) {
    if (dist_grid->vertex_owners[i].count == -1) {
 
      vertex_reorder[j] = i;
      ++j;
      coord_indices[offset] = i;
      cve_reorder[offset] = offset;
      ++offset;
    }
  }
  for (size_t i = 0, j = 0; i < dist_grid->num_edges; ++i) {
    if (dist_grid->edge_owners[i].count == -1) {
      // compute middle point of edge
      size_t * vertex_idxes = &(dist_grid->edge_to_vertex[i][0]);
      yac_int vertex_ids[2] = {dist_grid->vertex_ids[vertex_idxes[0]],
                               dist_grid->vertex_ids[vertex_idxes[1]]};
      edge_reorder[j] = i;
      ++j;
      coord_indices[offset] =
        vertex_idxes[(vertex_ids[0] > vertex_ids[1])?1:0];
      cve_reorder[offset] = offset;
      ++offset;
    }
  }
 
  yac_quicksort_index_size_t_size_t(
    coord_indices, total_num_missing, cve_reorder);
 
  size_t num_unique_idxs = 0;
  for (size_t i = 0, prev_idx = SIZE_MAX; i < total_num_missing; ++i) {
    size_t curr_idx = coord_indices[i];
    if (curr_idx != prev_idx) {
      prev_idx = curr_idx;
      ++num_unique_idxs;
    }
  }
 
  coordinate_pointer search_coords =
    xmalloc(num_unique_idxs * sizeof(*search_coords));
  int * search_ranks = xmalloc(num_unique_idxs * sizeof(*search_ranks));
 
  for (size_t i = 0, j = 0, prev_idx = SIZE_MAX; i < total_num_missing; ++i) {
    size_t curr_idx = coord_indices[i];
    if (curr_idx != prev_idx) {
      memcpy(search_coords[j], dist_grid->vertex_coordinates[curr_idx],
             3 * sizeof(double));
      prev_idx = curr_idx;
      ++j;
    }
    ranks_idxs[i] = j-1;
  }
 
  yac_proc_sphere_part_do_point_search(
    proc_sphere_part, search_coords, num_unique_idxs, search_ranks);
  free(search_coords);
 
  yac_quicksort_index_size_t_size_t(
    cve_reorder, total_num_missing, ranks_idxs);
 
  offset = 0;
  for (size_t i = 0; i < num_missing_vertices; ++i, ++offset)
      vertex_ranks[i] = search_ranks[ranks_idxs[offset]];
  for (size_t i = 0; i < num_missing_edges; ++i, ++offset)
      edge_ranks[i] = search_ranks[ranks_idxs[offset]];
  free(search_ranks);
 
  yac_quicksort_index_int_size_t(
    vertex_ranks, num_missing_vertices, vertex_reorder);
  yac_quicksort_index_int_size_t(
    edge_ranks, num_missing_edges, edge_reorder);
 
  int * sendcounts, * recvcounts, * sdispls, * rdispls;
  yac_get_comm_buffers(
    2, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
  int * int_buffer =
    xmalloc(4 * (size_t)comm_size * sizeof(*int_buffer));
  int * total_sendcounts = int_buffer + 0 * comm_size;
  int * total_recvcounts = int_buffer + 1 * comm_size;
  int * total_sdispls =    int_buffer + 2 * comm_size;
  int * total_rdispls =    int_buffer + 3 * comm_size;
 
  for (size_t i = 0; i < num_missing_vertices; ++i)
      sendcounts[2 * vertex_ranks[i] + 0]++;
  for (size_t i = 0; i < num_missing_edges; ++i)
      sendcounts[2 * edge_ranks[i] + 1]++;
 
  yac_generate_alltoallv_args(
    2, sendcounts, recvcounts, sdispls, rdispls, comm);
 
  size_t num_requested_ids[2] = {0,0};
  for (int i = 0, saccu = 0, raccu = 0; i < comm_size; ++i) {
    total_sdispls[i] = saccu;
    total_rdispls[i] = raccu;
    total_sendcounts[i] = 0;
    total_recvcounts[i] = 0;
    for (int j = 0; j < 2; ++j) {
      total_sendcounts[i] += sendcounts[2 * i + j];
      total_recvcounts[i] += recvcounts[2 * i + j];
      num_requested_ids[j] += (size_t)(recvcounts[2 * i + j]);
    }
    saccu += total_sendcounts[i];
    raccu += total_recvcounts[i];
  }
  size_t request_count = (size_t)(total_recvcounts[comm_size - 1]) +
                         (size_t)(total_rdispls[comm_size - 1]);
 
  yac_int * exchange_id_buffer =
    xmalloc(
      (total_num_missing + 2 * request_count) * sizeof(*exchange_id_buffer));
  yac_int * id_send_buffer = exchange_id_buffer + request_count;
  yac_int * id_recv_buffer = exchange_id_buffer + request_count + total_num_missing;
 
  // pack the ids for which the local process is missing the owner information
  for (size_t i = 0; i < num_missing_vertices; ++i) {
    int pos = sdispls[2 * vertex_ranks[i] + 0 + 1]++;
    id_send_buffer[pos] = dist_grid->vertex_ids[vertex_reorder[i]];
  }
  for (size_t i = 0; i < num_missing_edges; ++i) {
    int pos = sdispls[2 * edge_ranks[i] + 1 + 1]++;
    id_send_buffer[pos] = dist_grid->edge_ids[edge_reorder[i]];
  }
 
  // exchange the ids for which the local process is missing the owner information
  yac_alltoallv_yac_int_p2p(
    id_send_buffer, total_sendcounts, total_sdispls,
    id_recv_buffer, total_recvcounts, total_rdispls, comm);
 
  yac_int * requested_ids[2] =
    {exchange_id_buffer, exchange_id_buffer + num_requested_ids[0]};
 
  // unpack data
  {
    size_t ids_offset[2] = {0, 0}, offset = 0;
    for (int i = 0; i < comm_size; ++i) {
      for (int j = 0; j < 2; ++j) {
        size_t curr_count = recvcounts[2 * i + j];
        memcpy(requested_ids[j] + ids_offset[j], id_recv_buffer + offset,
               curr_count * sizeof(*(requested_ids[0])));
        ids_offset[j] += curr_count;
        offset += curr_count;
      }
    }
  }
 
  struct remote_point_infos ** remote_point_infos_buffer =
    xmalloc(request_count * sizeof(*remote_point_infos_buffer));
  struct remote_point_infos ** found_points[2] =
    {remote_point_infos_buffer,
     remote_point_infos_buffer + num_requested_ids[0]};
  struct remote_points * local_remote_points[2] =
    {dist_vertex_owner, dist_edge_owner};
 
  size_t max_num_requested_ids = MAX(num_requested_ids[0], num_requested_ids[1]);
  size_t * reorder_buffer = xmalloc(max_num_requested_ids * sizeof(*reorder_buffer));
  // lookup owners for requested ids
  for (int i = 0; i < 2; ++i)
    lookup_remote_points(
      requested_ids[i], found_points[i], reorder_buffer, num_requested_ids[i],
      local_remote_points[i]);
  free(reorder_buffer);
  free(exchange_id_buffer);
 
  int * pack_size_buffer =
    xmalloc((total_num_missing + request_count) * sizeof(*pack_size_buffer));
  int * send_pack_sizes = pack_size_buffer;
  int * recv_pack_sizes = pack_size_buffer + request_count;
  MPI_Datatype point_info_dt = yac_get_remote_point_info_mpi_datatype(comm);
 
  {
    size_t offset = 0, offsets[2] = {0, 0};
    for (int i = 0; i < comm_size; ++i) {
      for (int j = 0; j < 2; ++j) {
        size_t curr_count = (size_t)(recvcounts[2 * i + j]);
        // get pack sizes
        yac_remote_point_infos_get_pack_sizes(
          found_points[j] + offsets[j], curr_count,
          send_pack_sizes + offset, point_info_dt, comm);
        offsets[j] += curr_count;
        offset += curr_count;
      }
    }
  }
 
  // exchange pack sizes
  yac_alltoallv_int_p2p(
    send_pack_sizes, total_recvcounts, total_rdispls,
    recv_pack_sizes, total_sendcounts, total_sdispls, comm);
 
  size_t send_buffer_size = 0, recv_buffer_size = 0;
  for (size_t i = 0; i < request_count; ++i)
    send_buffer_size += (size_t)(send_pack_sizes[i]);
  for (size_t i = 0; i < total_num_missing; ++i)
    recv_buffer_size += (size_t)(recv_pack_sizes[i]);
 
  void * result_buffer = xmalloc(send_buffer_size + recv_buffer_size);
  void * send_result_buffer = result_buffer;
  void * recv_result_buffer = (unsigned char *)result_buffer + send_buffer_size;
 
  // pack data
  {
    size_t send_buffer_size = 0, offset = 0, offsets[2] = {0,0};
    for (int i = 0; i < comm_size; ++i) {
      for (int j = 0; j < 2; ++j) {
        size_t curr_count = (size_t)(recvcounts[2 * i + j]);
        // get pack sizes
        yac_remote_point_infos_pack(
          found_points[j] + offsets[j],
          send_pack_sizes + offset, curr_count,
          (void*)((unsigned char*)send_result_buffer + send_buffer_size),
          point_info_dt, comm);
        for (size_t k = 0; k < curr_count; ++k)
          send_buffer_size += send_pack_sizes[offset + k];
        offsets[j] += curr_count;
        offset += curr_count;
      }
    }
  }
 
  {
    size_t send_offset = 0, recv_offset = 0;
    size_t saccu = 0, raccu = 0;
    for (int i = 0; i < comm_size; ++i) {
      total_sdispls[i] = (int)saccu;
      total_rdispls[i] = (int)raccu;
      int curr_sendcount = 0, curr_recvcount = 0;
      for (int k = 0; k < 2; ++k) {
        for (int j = 0; j < recvcounts[2 * i + k]; ++j, ++send_offset)
          curr_sendcount += (size_t)(send_pack_sizes[send_offset]);
        for (int j = 0; j < sendcounts[2 * i + k]; ++j, ++recv_offset)
          curr_recvcount += (size_t)(recv_pack_sizes[recv_offset]);
      }
      total_sendcounts[i] = (int)curr_sendcount;
      total_recvcounts[i] = (int)curr_recvcount;
 
      YAC_ASSERT(
        (saccu <= INT_MAX) && (raccu <= INT_MAX) &&
        (curr_sendcount <= INT_MAX) && (curr_recvcount <= INT_MAX),
        "ERROR(generate_dist_grid_remote_point_infos): "
        "int argument to MPI is bigger than INT_MAX")
 
      saccu += curr_sendcount;
      raccu += curr_recvcount;
    }
  }
 
  // exchange results
  yac_alltoallv_p2p(
    send_result_buffer, total_sendcounts, total_sdispls,
    recv_result_buffer, total_recvcounts, total_rdispls, 1, MPI_PACKED, comm);
 
  // unpack results
  {
    size_t counts[2] = {0, 0}, count = 0, offset = 0;
    struct remote_point_infos * remote_points[2] =
      {dist_grid->vertex_owners, dist_grid->edge_owners};
    size_t * reorder[2] = {vertex_reorder, edge_reorder};
    for (int i = 0; i < comm_size; ++i) {
      for (int j = 0; j < 2; ++j) {
        size_t curr_count = sendcounts[2 * i + j];
        for (size_t k = 0; k < curr_count; ++k, ++counts[j], ++count) {
          int curr_buffer_size =
            recv_pack_sizes[(size_t)(sdispls[2 * i + j]) + k];
          int position = 0;
          yac_remote_point_infos_unpack(
            (void*)((unsigned char*)recv_result_buffer + offset),
            curr_buffer_size, &position,
            remote_points[j] + reorder[j][counts[j]], point_info_dt, comm);
          offset += (size_t)(curr_buffer_size);
        }
      }
    }
  }
 
  yac_mpi_call(MPI_Type_free(&point_info_dt), comm);
 
  yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
  free(int_buffer);
  free(result_buffer);
  free(pack_size_buffer);
  free(remote_point_infos_buffer);
  free(size_t_buffer);
  free(ranks_buffer);
}
 
static int compare_remote_point(const void * a, const void * b) {
 
  return ((const struct remote_point*)a)->global_id -
         ((const struct remote_point*)b)->global_id;
}
 
// inserts an element into an array and increases the corresponding size
// if the element already exists in the array nothing is done
static void insert_global_id(yac_int * ids, size_t n, yac_int id) {
 
  size_t i;
  for (i = 0; i < n; ++i) if (ids[i] >= id) break;
  // copy new id into array and move bigger elements one position up
  if (n != i) memmove(ids + i + 1, ids + i, (n - i) * sizeof(*ids));
  ids[i] = id;
}
 
// inserts an element into an array and increases the corresponding size
// if the element already exists in the array nothing is done
static void insert_rank(int * ranks, int * count, int rank) {
 
  int i;
  int n = *count;
 
  for (i = 0; i < n; ++i) if (ranks[i] >= rank) break;
 
  // if the rank is already in the array
  if (i != n) {
    if (ranks[i] == rank) return;
    else memmove(ranks + i + 1, ranks + i, ((size_t)(n - i)) * sizeof(*ranks));
  }
  ranks[i] = rank;
  *count = n + 1;
}
 
static int compare_n_ids_reorder_ids (
  const void * a, const void * b) {
 
  int count_a = ((struct n_ids_reorder const *)a)->count;
  int count_b = ((struct n_ids_reorder const *)b)->count;
  yac_int * a_ids = ((struct n_ids_reorder const *)a)->ids;
  yac_int * b_ids = ((struct n_ids_reorder const *)b)->ids;
  int ret = count_a - count_b;
  for (int i = 0; !ret && (i < count_a); ++i)
    ret = (a_ids[i] > b_ids[i]) - (a_ids[i] < b_ids[i]);
  return ret;
}
 
static int compare_n_ids_reorder_reorder (
  const void * a, const void * b) {
 
  return (((struct n_ids_reorder const *)a)->reorder_idx >
          ((struct n_ids_reorder const *)b)->reorder_idx) -
         (((struct n_ids_reorder const *)a)->reorder_idx <
          ((struct n_ids_reorder const *)b)->reorder_idx);
}
 
static void determine_dist_cell_ranks(
  struct proc_sphere_part_node * proc_sphere_part,
  struct basic_grid_data * grid, MPI_Comm comm, int ** dist_cell_ranks,
  int ** dist_cell_rank_counts, size_t ** dist_cell_rank_offsets) {
 
  int comm_size;
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  size_t num_cells = grid->num_cells;
  int * ranks_buffer = xmalloc((size_t)comm_size * sizeof(*ranks_buffer));
  size_t dist_cell_ranks_array_size = num_cells;
  int * dist_cell_ranks_ = xmalloc(num_cells * sizeof(*dist_cell_ranks_));
  int * dist_cell_rank_counts_ =
    (*dist_cell_rank_counts =
       xmalloc(num_cells * sizeof(*dist_cell_rank_counts_)));
  size_t * dist_cell_rank_offsets_ =
    (*dist_cell_rank_offsets =
       xmalloc(num_cells * sizeof(*dist_cell_rank_offsets_)));
  size_t offset = 0;
 
  // determine maximum  number of vertices per cell
  int max_num_vertices_per_cell = 0;
  for (size_t i = 0; i < num_cells; ++i)
    if (grid->num_vertices_per_cell[i] > max_num_vertices_per_cell)
      max_num_vertices_per_cell = grid->num_vertices_per_cell[i];
 
  struct grid_cell cell;
  cell.coordinates_xyz =
    xmalloc(
      (size_t)max_num_vertices_per_cell * sizeof(*(cell.coordinates_xyz)));
  cell.edge_type =
    xmalloc((size_t)max_num_vertices_per_cell * sizeof(*(cell.edge_type)));
 
  size_t * cell_to_vertex = grid->cell_to_vertex;
  size_t * cell_to_vertex_offsets = grid->cell_to_vertex_offsets;
  size_t * cell_to_edge = grid->cell_to_edge;
  size_t * cell_to_edge_offsets = grid->cell_to_edge_offsets;
  coordinate_pointer vertex_coordinates = grid->vertex_coordinates;
  enum yac_edge_type * edge_type = grid->edge_type;
  int * num_vertices_per_cell = grid->num_vertices_per_cell;
 
  for (size_t i = 0; i < num_cells; ++i) {
 
    // get bounding circle for current cell
    cell.num_corners = num_vertices_per_cell[i];
    size_t * curr_cell_to_vertex =
      cell_to_vertex + cell_to_vertex_offsets[i];
    size_t * curr_cell_to_edge =
      cell_to_edge + cell_to_edge_offsets[i];
    for (int j = 0; j < num_vertices_per_cell[i]; ++j) {
      coordinate_pointer curr_vertex_coords =
        vertex_coordinates + curr_cell_to_vertex[j];
      for (int k = 0; k < 3; ++k)
        cell.coordinates_xyz[j][k] = (*curr_vertex_coords)[k];
      cell.edge_type[j] = edge_type[curr_cell_to_edge[j]];
    }
    int rank_count;
    if (cell.num_corners > 0) {
 
      struct bounding_circle bnd_circle;
      yac_get_cell_bounding_circle(cell, &bnd_circle);
 
      yac_proc_sphere_part_do_bnd_circle_search(
        proc_sphere_part, bnd_circle, ranks_buffer, &rank_count);
    } else {
      ranks_buffer[0] = 0;
      rank_count = 1;
    }
 
    ENSURE_ARRAY_SIZE(dist_cell_ranks_, dist_cell_ranks_array_size,
                      offset + (size_t)rank_count);
    memcpy(dist_cell_ranks_ + offset, ranks_buffer,
           (size_t)rank_count * sizeof(*ranks_buffer));
    dist_cell_rank_counts_[i] = rank_count;
    dist_cell_rank_offsets_[i] = offset;
    offset += (size_t)rank_count;
  }
 
  free(cell.edge_type);
  free(cell.coordinates_xyz);
  free(ranks_buffer);
 
  *dist_cell_ranks =
    xrealloc(dist_cell_ranks_, offset * sizeof(*dist_cell_ranks_));
}
 
static MPI_Datatype yac_get_id_reorder_coord_coord_mpi_datatype(MPI_Comm comm) {
 
  struct id_reorder_coord dummy;
  MPI_Datatype coord_dt;
  int array_of_blocklengths[] = {3};
  const MPI_Aint array_of_displacements[] =
    {(MPI_Aint)((void*)&(dummy.coord) - (void*)&dummy)};
  const MPI_Datatype array_of_types[] = {MPI_DOUBLE};
  yac_mpi_call(
    MPI_Type_create_struct(1, array_of_blocklengths, array_of_displacements,
                           array_of_types, &coord_dt), comm);
  return yac_create_resized(coord_dt, sizeof(dummy), comm);
}
 
static MPI_Datatype yac_get_id_reorder_coord_id_mpi_datatype(MPI_Comm comm) {
 
  struct id_reorder_coord dummy;
  MPI_Datatype global_id_dt;
  int array_of_blocklengths[] = {1};
  const MPI_Aint array_of_displacements[] =
    {(MPI_Aint)((void*)&(dummy.global_id) - (void*)&dummy)};
  const MPI_Datatype array_of_types[] = {Xt_int_dt};
  yac_mpi_call(
    MPI_Type_create_struct(1, array_of_blocklengths, array_of_displacements,
                           array_of_types, &global_id_dt), comm);
  return yac_create_resized(global_id_dt, sizeof(dummy), comm);
}
 
static int compare_id_reorder_coord_coord (
  const void * a, const void * b) {
 
  return compare_coords(((const struct id_reorder_coord *)a)->coord,
                        ((const struct id_reorder_coord *)b)->coord);
}
 
static int compare_id_reorder_coord_reorder_idx (
  const void * a, const void * b) {
 
  return (((const struct id_reorder_coord *)a)->reorder_idx >
          ((const struct id_reorder_coord *)b)->reorder_idx) -
         (((const struct id_reorder_coord *)a)->reorder_idx <
          ((const struct id_reorder_coord *)b)->reorder_idx);
}
 
static int * generate_vertex_ids(
  struct proc_sphere_part_node * proc_sphere_part,
  struct basic_grid_data * grid, MPI_Comm comm) {
 
  int * vertex_ranks = xmalloc(grid->num_vertices * sizeof(*vertex_ranks));
 
  // determine the dist grid ranks for all vertices
  yac_proc_sphere_part_do_point_search(
    proc_sphere_part, grid->vertex_coordinates, grid->num_vertices,
    vertex_ranks);
 
  // check whether only of subset of the processes have defined
  // their global ids, which is not supported
  int ids_available_local =
    (grid->num_vertices > 0) && (grid->vertex_ids != NULL);
  int ids_available_global;
  yac_mpi_call(MPI_Allreduce(&ids_available_local, &ids_available_global, 1,
                             MPI_INT, MPI_MAX, comm), comm);
 
  // if there are vertices defined locally for the current grid
  if (grid->num_vertices != 0) {
    YAC_ASSERT(
      ids_available_local || !ids_available_global,
      "ERROR(generate_dist_remote_points_vertex): inconsistent global ids")
 
    // generate core mask if not available
    if (grid->core_vertex_mask == NULL) {
      int * core_mask =
        (grid->core_vertex_mask =
           xmalloc(grid->num_vertices * sizeof(*core_mask)));
      for (size_t j = 0; j < grid->num_vertices; ++j) core_mask[j] = 1;
    }
  }
 
  // if we do not need to generate the global ids
  if (ids_available_global) return vertex_ranks;
 
  int comm_size;
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  int * sendcounts, * recvcounts, * sdispls, * rdispls;
  yac_get_comm_buffers(
    1, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
 
  for (size_t i = 0; i < grid->num_vertices; ++i)
    sendcounts[vertex_ranks[i]]++;
 
  yac_generate_alltoallv_args(
    1, sendcounts, recvcounts, sdispls, rdispls, comm);
 
  size_t orig_vertex_count = grid->num_vertices;
  size_t dist_vertex_count = (size_t)(recvcounts[comm_size - 1]) +
                             (size_t)(rdispls[comm_size - 1]);
 
  struct id_reorder_coord * id_reorder_coord_buffer =
    xmalloc((orig_vertex_count + dist_vertex_count) *
            sizeof(*id_reorder_coord_buffer));
  struct id_reorder_coord *
    id_reorder_coord_send_buffer = id_reorder_coord_buffer;
  struct id_reorder_coord *
    id_reorder_coord_recv_buffer = id_reorder_coord_buffer + orig_vertex_count;
 
  // pack send buffer
  for (size_t i = 0; i < orig_vertex_count; ++i) {
    int pos = sdispls[vertex_ranks[i] + 1]++;
    id_reorder_coord_send_buffer[pos].reorder_idx = i;
    for (int j = 0; j < 3; ++j)
      id_reorder_coord_send_buffer[pos].coord[j] =
        grid->vertex_coordinates[i][j];
  }
 
  MPI_Datatype id_reorder_coord_coord_dt =
    yac_get_id_reorder_coord_coord_mpi_datatype(comm);
 
  // exchange data
  yac_alltoallv_p2p(
    id_reorder_coord_send_buffer, sendcounts, sdispls,
    id_reorder_coord_recv_buffer, recvcounts, rdispls,
    sizeof(*id_reorder_coord_send_buffer), id_reorder_coord_coord_dt, comm);
 
  yac_mpi_call(MPI_Type_free(&id_reorder_coord_coord_dt), comm);
 
  for (size_t i = 0; i < dist_vertex_count; ++i)
    id_reorder_coord_recv_buffer[i].reorder_idx = i;
 
  // sort received vertices based on coordinates
  qsort(id_reorder_coord_recv_buffer, dist_vertex_count,
        sizeof(*id_reorder_coord_recv_buffer), compare_id_reorder_coord_coord);
 
  size_t unique_count = dist_vertex_count > 0;
  struct id_reorder_coord * prev = id_reorder_coord_recv_buffer;
 
  // determine number of unique coordinates
  for (size_t i = 0; i < dist_vertex_count; ++i) {
    struct id_reorder_coord * curr = id_reorder_coord_recv_buffer + i;
    if (compare_id_reorder_coord_coord(prev, curr)) {
      ++unique_count;
      prev = curr;
    }
    curr->global_id = (yac_int)unique_count - 1;
  }
 
  YAC_ASSERT(
    unique_count <= (size_t)XT_INT_MAX,
    "ERROR(generate_vertex_ids): global_id out of bounds")
 
  yac_int yac_int_unique_count = (yac_int)unique_count;
  yac_int id_offset;
 
  // determine exclusive scan of sum of numbers of unique
  // coordinates on all ranks
  yac_mpi_call(MPI_Exscan(&yac_int_unique_count, &id_offset, 1, Xt_int_dt,
                          MPI_SUM, comm), comm);
  int comm_rank;
  yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
  if (comm_rank == 0) id_offset = 0;
 
  YAC_ASSERT(
    ((size_t)id_offset + unique_count) <= (size_t)XT_INT_MAX,
    "ERROR(generate_vertex_ids): global_id out of bounds")
 
  // adjust global ids
  for (size_t i = 0; i < dist_vertex_count; ++i)
    id_reorder_coord_recv_buffer[i].global_id += id_offset;
 
  // return received vertices into original order
  qsort(id_reorder_coord_recv_buffer, dist_vertex_count,
        sizeof(*id_reorder_coord_recv_buffer),
        compare_id_reorder_coord_reorder_idx);
 
  MPI_Datatype id_reorder_coord_id_dt =
    yac_get_id_reorder_coord_id_mpi_datatype(comm);
 
  // return generated global ids data
  yac_alltoallv_p2p(
    id_reorder_coord_recv_buffer, recvcounts, rdispls,
    id_reorder_coord_send_buffer, sendcounts, sdispls,
    sizeof(*id_reorder_coord_send_buffer), id_reorder_coord_id_dt, comm);
 
  yac_mpi_call(MPI_Type_free(&id_reorder_coord_id_dt), comm);
 
  yac_int * vertex_ids =
    ((grid->vertex_ids =
      (grid->num_vertices > 0)?
        xmalloc(grid->num_vertices * sizeof(*(grid->vertex_ids))):NULL));
 
  for (size_t i = 0; i < grid->num_vertices; ++i)
    vertex_ids[id_reorder_coord_send_buffer[i].reorder_idx] =
      id_reorder_coord_send_buffer[i].global_id;
 
  free(id_reorder_coord_buffer);
  yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
 
  return vertex_ranks;
}
 
static void generate_ce_ids(
  struct proc_sphere_part_node * proc_sphere_part,
  struct basic_grid_data * grid, int * vertex_ranks, MPI_Comm comm) {
 
  int comm_rank, comm_size;
  yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  // check whether only of subset of the processes have defined
  // their global ids, which is not supported
  int ids_available_local[2], ids_available_global[2];
  ids_available_local[0] = (grid->num_cells > 0) && (grid->cell_ids != NULL);
  ids_available_local[1] = (grid->num_edges > 0) && (grid->edge_ids != NULL);
  yac_mpi_call(MPI_Allreduce(ids_available_local, ids_available_global, 2,
                             MPI_INT, MPI_MAX, comm), comm);
 
  if (grid->num_cells > 0) {
    YAC_ASSERT(
      ids_available_local[0] == ids_available_global[0],
      "ERROR(generate_ce_ids): inconsistent global ids")
    // generate core mask if not available
    if (grid->core_cell_mask == NULL) {
      int * core_cell_mask =
        (grid->core_cell_mask =
           xmalloc(grid->num_cells * sizeof(*core_cell_mask)));
      for (size_t j = 0; j < grid->num_cells; ++j) core_cell_mask[j] = 1;
    }
  }
  if (grid->num_edges > 0) {
    YAC_ASSERT(
      ids_available_local[1] == ids_available_global[1],
      "ERROR(generate_ce_ids): inconsistent global ids")
    // generate core mask if not available
    if (grid->core_edge_mask == NULL) {
      int * core_edge_mask =
        (grid->core_edge_mask =
           xmalloc(grid->num_edges * sizeof(*core_edge_mask)));
      for (size_t j = 0; j < grid->num_edges; ++j) core_edge_mask[j] = 1;
    }
  }
 
  // if no ids need to be generated
  if (ids_available_global[0] && ids_available_global[1]) return;
 
  int * rank_buffer =
    xmalloc(
      (((ids_available_global[0])?0:(grid->num_cells)) +
       ((ids_available_global[1])?0:(grid->num_edges))) * sizeof(*rank_buffer));
  int * cell_ranks = rank_buffer;
  int * edge_ranks =
    rank_buffer + ((ids_available_global[0])?0:(grid->num_cells));
 
  int * int_buffer =
    xmalloc((8 * (size_t)comm_size + 1) * sizeof(*int_buffer));
  int * sendcounts = int_buffer + 0 * comm_size;
  int * recvcounts = int_buffer + 2 * comm_size;
  int * total_sendcounts = int_buffer + 4 * comm_size;
  int * total_recvcounts = int_buffer + 5 * comm_size;
  int * total_sdispls = int_buffer + 6 * comm_size;
  int * total_rdispls = int_buffer + 7 * comm_size + 1;
  memset(sendcounts, 0, 2 * (size_t)comm_size * sizeof(*sendcounts));
 
  yac_int * cell_to_vertex_ids = NULL;
  int max_num_vertices_per_cell = 0;
 
  if (!ids_available_global[0]) {
 
    { // compute the global maximum number of vertices per cell
      for (size_t i = 0; i < grid->num_cells; ++i)
        if (grid->num_vertices_per_cell[i] > max_num_vertices_per_cell)
          max_num_vertices_per_cell = grid->num_vertices_per_cell[i];
      yac_mpi_call(MPI_Allreduce(MPI_IN_PLACE, &max_num_vertices_per_cell, 1,
                                 MPI_INT, MPI_MAX, comm), comm);
    }
 
    cell_to_vertex_ids =
      xmalloc(grid->num_cells * max_num_vertices_per_cell *
              sizeof(*cell_to_vertex_ids));
 
    for (size_t i = 0; i < grid->num_cells; ++i) {
 
      size_t curr_num_vertices = grid->num_vertices_per_cell[i];
      yac_int * curr_cell_to_vertex_ids =
        cell_to_vertex_ids + i * max_num_vertices_per_cell;
 
      int cell_rank;
      if (curr_num_vertices > 0) {
        size_t * curr_cell_vertices =
          grid->cell_to_vertex + grid->cell_to_vertex_offsets[i];
        size_t min_vertex = curr_cell_vertices[0];
        curr_cell_to_vertex_ids[0] = grid->vertex_ids[min_vertex];
        for (size_t j = 1; j < curr_num_vertices; ++j) {
          size_t curr_vertex_idx = curr_cell_vertices[j];
          yac_int curr_vertex_id = grid->vertex_ids[curr_vertex_idx];
          insert_global_id(curr_cell_to_vertex_ids, j, curr_vertex_id);
          if (curr_cell_to_vertex_ids[0] == curr_vertex_id)
            min_vertex = curr_vertex_idx;
        }
        cell_rank = vertex_ranks[min_vertex];
      } else {
        cell_rank = 0;
      }
      for (size_t j = curr_num_vertices; j < max_num_vertices_per_cell; ++j)
        curr_cell_to_vertex_ids[j] = XT_INT_MAX;
 
      sendcounts[2 * ((cell_ranks[i] = cell_rank)) + 0]++;
    }
  }
 
  yac_int * edge_to_vertex_ids = NULL;
 
  if (!ids_available_global[1]) {
 
    edge_to_vertex_ids =
      xmalloc(2 * grid->num_edges * sizeof(*edge_to_vertex_ids));
 
    for (size_t i = 0; i < grid->num_edges; ++i) {
 
      size_t * curr_edge_to_vertex = &(grid->edge_to_vertex[i][0]);
      yac_int * curr_edge_vertex_ids = edge_to_vertex_ids + 2 * i;
      curr_edge_vertex_ids[0] = grid->vertex_ids[curr_edge_to_vertex[0]];
      curr_edge_vertex_ids[1] = grid->vertex_ids[curr_edge_to_vertex[1]];
 
      if (curr_edge_vertex_ids[0] > curr_edge_vertex_ids[1]) {
        yac_int temp = curr_edge_vertex_ids[0];
        curr_edge_vertex_ids[0] = curr_edge_vertex_ids[1];
        curr_edge_vertex_ids[1] = temp;
        sendcounts[
          2 * ((edge_ranks[i] = vertex_ranks[curr_edge_to_vertex[1]])) + 1]++;
      } else {
        sendcounts[
          2 * ((edge_ranks[i] = vertex_ranks[curr_edge_to_vertex[0]])) + 1]++;
      }
    }
  }
 
  // exchange the number of cells and edges
  yac_mpi_call(MPI_Alltoall(sendcounts, 2, MPI_INT,
                            recvcounts, 2, MPI_INT, comm), comm);
 
  total_sdispls[0] = 0;
  size_t recv_counts[2] = {0,0};
  for (int i = 0, saccu = 0, raccu = 0; i < comm_size; ++i) {
    total_sdispls[i+1] = saccu;
    total_rdispls[i] = raccu;
    recv_counts[0] += recvcounts[2 * i + 0];
    recv_counts[1] += recvcounts[2 * i + 1];
    total_sendcounts[i] = sendcounts[2 * i + 0] * max_num_vertices_per_cell +
                          sendcounts[2 * i + 1] * 2;
    total_recvcounts[i] = recvcounts[2 * i + 0] * max_num_vertices_per_cell +
                          recvcounts[2 * i + 1] * 2;
    saccu += total_sendcounts[i];
    raccu += total_recvcounts[i];
  }
  size_t local_data_count = (size_t)(total_sendcounts[comm_size - 1]) +
                            (size_t)(total_sdispls[comm_size]);
  size_t recv_count = (size_t)(total_recvcounts[comm_size - 1]) +
                      (size_t)(total_rdispls[comm_size - 1]);
 
  yac_int * yac_int_buffer =
    xmalloc((local_data_count + recv_count) * sizeof(*yac_int_buffer));
  yac_int * send_buffer = yac_int_buffer;
  yac_int * recv_buffer = yac_int_buffer + local_data_count;
 
  // pack send buffer
  if (!ids_available_global[0])
    for (size_t i = 0; i < grid->num_cells; ++i)
      for (int j = 0; j < max_num_vertices_per_cell; ++j)
        send_buffer[total_sdispls[cell_ranks[i] + 1]++] =
          cell_to_vertex_ids[i * max_num_vertices_per_cell + j];
  if (!ids_available_global[1])
    for (size_t i = 0; i < grid->num_edges; ++i)
      for (int j = 0; j < 2; ++j)
        send_buffer[total_sdispls[edge_ranks[i] + 1]++] =
          edge_to_vertex_ids[2 * i + j];
 
  free(edge_to_vertex_ids);
  free(cell_to_vertex_ids);
 
  // exchange data
  yac_alltoallv_p2p(
    send_buffer, total_sendcounts, total_sdispls,
    recv_buffer, total_recvcounts, total_rdispls,
    sizeof(*send_buffer), Xt_int_dt, comm);
 
  struct n_ids_reorder * n_ids_reorder_buffer =
    xmalloc((recv_counts[0] + recv_counts[1]) * sizeof(*n_ids_reorder_buffer));
  struct n_ids_reorder * n_ids_reorder[2] =
    {n_ids_reorder_buffer, n_ids_reorder_buffer + recv_counts[0]};
 
  size_t offset = 0;
  int index_counts[2] = {max_num_vertices_per_cell, 2};
  size_t reorder_idx = 0;
  recv_counts[0] = 0;
  recv_counts[1] = 0;
  for (int i = 0; i < comm_size; ++i) {
    for (int j = 0; j < 2; ++j) {
      int curr_count = recvcounts[2 * i + j];
      for (int k = 0; k < curr_count; ++k, ++reorder_idx, ++recv_counts[j]) {
        n_ids_reorder[j][recv_counts[j]].count = index_counts[j];
        n_ids_reorder[j][recv_counts[j]].ids = recv_buffer + offset;
        n_ids_reorder[j][recv_counts[j]].reorder_idx = reorder_idx;
        offset += index_counts[j];
      }
    }
  }
 
  for (int i = 0; i < 2; ++i) {
 
    if (ids_available_global[i]) continue;
 
    qsort(n_ids_reorder[i], recv_counts[i], sizeof(*(n_ids_reorder[i])),
          compare_n_ids_reorder_ids);
 
    size_t unique_count = recv_counts[i] > 0;
    struct n_ids_reorder * prev = n_ids_reorder[i];
    struct n_ids_reorder * curr = n_ids_reorder[i];
 
    for (size_t j = 0; j < recv_counts[i]; ++j, ++curr) {
      if (compare_n_ids_reorder_ids(prev, curr)) {
        ++unique_count;
        prev = curr;
      }
      curr->global_id = (yac_int)(unique_count - 1);
    }
 
    YAC_ASSERT(
      unique_count <= (size_t)XT_INT_MAX,
      "ERROR(generate_global_ce_ids): global_id out of bounds")
 
    yac_int yac_int_unique_count = (yac_int)unique_count;
    yac_int id_offset;
 
    // determine exclusive scan of sum of numbers of unique ids on all ranks
    yac_mpi_call(MPI_Exscan(&yac_int_unique_count, &id_offset, 1, Xt_int_dt,
                            MPI_SUM, comm), comm);
    if (comm_rank == 0) id_offset = 0;
 
    YAC_ASSERT(
      ((size_t)id_offset + unique_count) <= (size_t)XT_INT_MAX,
      "ERROR(generate_global_ce_ids): global_id out of bounds")
 
    // adjust global ids
    for (size_t j = 0; j < recv_counts[i]; ++j)
      n_ids_reorder[i][j].global_id += id_offset;
  }
  free(yac_int_buffer);
 
  qsort(n_ids_reorder_buffer, recv_counts[0] + recv_counts[1],
        sizeof(*n_ids_reorder_buffer), compare_n_ids_reorder_reorder);
 
  yac_int * global_ids_buffer =
    xmalloc((recv_counts[0] + recv_counts[1] +
             ((ids_available_global[0])?0:(grid->num_cells)) +
             ((ids_available_global[1])?0:(grid->num_edges))) *
             sizeof(*global_ids_buffer));
  yac_int * send_global_ids = global_ids_buffer;
  yac_int * recv_global_ids =
    global_ids_buffer + recv_counts[0] + recv_counts[1];
 
  for (size_t i = 0; i < recv_counts[0] + recv_counts[1]; ++i)
    send_global_ids[i] = n_ids_reorder_buffer[i].global_id;
  free(n_ids_reorder_buffer);
 
  // generate count and displs data
  for (int i = 0, saccu = 0, raccu = 0; i < comm_size; ++i) {
    total_sdispls[i] = saccu;
    total_rdispls[i] = raccu;
    saccu +=
      ((total_sendcounts[i] = recvcounts[2 * i + 0] + recvcounts[2 * i + 1]));
    raccu +=
      ((total_recvcounts[i] = sendcounts[2 * i + 0] + sendcounts[2 * i + 1]));
  }
 
  // exchange generated global ids data
  yac_alltoallv_p2p(
    send_global_ids, total_sendcounts, total_sdispls,
    recv_global_ids, total_recvcounts, total_rdispls,
    sizeof(*send_global_ids), Xt_int_dt, comm);
 
  if ((!ids_available_global[0]) && (grid->num_cells > 0))
    grid->cell_ids = xmalloc(grid->num_cells * sizeof(*grid->cell_ids));
  if ((!ids_available_global[1]) && (grid->num_edges > 0))
    grid->edge_ids = xmalloc(grid->num_edges * sizeof(grid->edge_ids));
 
  // unpack generated global ids
  if (!ids_available_global[0])
    for (size_t i = 0; i < grid->num_cells; ++i)
      grid->cell_ids[i] = recv_global_ids[total_rdispls[cell_ranks[i]]++];
  if (!ids_available_global[1])
    for (size_t i = 0; i < grid->num_edges; ++i)
      grid->edge_ids[i] = recv_global_ids[total_rdispls[edge_ranks[i]]++];
 
  free(rank_buffer);
  free(int_buffer);
  free(global_ids_buffer);
}
 
static void generate_global_ids(
  struct proc_sphere_part_node * proc_sphere_part,
  struct basic_grid_data * grid, MPI_Comm comm) {
 
  // generate global ids for all vertices (if non-existent)
  int * vertex_ranks = generate_vertex_ids(proc_sphere_part, grid,  comm);
 
  // generate global ids for all cell and edge (if non-existent)
  generate_ce_ids(proc_sphere_part, grid, vertex_ranks, comm);
 
  free(vertex_ranks);
}
 
static size_t_2_pointer generate_edge_to_cell(
  const_size_t_pointer cell_to_edge, const_int_pointer num_edges_per_cell,
  int * core_cell_mask, size_t num_cells, size_t num_edges) {
 
  if (num_cells == 0) return NULL;
 
  size_t_2_pointer edge_to_cell = xmalloc(num_edges * sizeof(*edge_to_cell));
 
  for (size_t i = 0; i < num_edges; ++i) {
    edge_to_cell[i][0] = SIZE_MAX;
    edge_to_cell[i][1] = SIZE_MAX;
  }
 
  for (size_t i = 0, offset = 0; i < num_cells; ++i) {
 
    size_t curr_num_edges = num_edges_per_cell[i];
    const_size_t_pointer curr_cell_to_edge = cell_to_edge + offset;
    offset += curr_num_edges;
 
    if ((core_cell_mask == NULL) || core_cell_mask[i]) {
 
      for (size_t j = 0; j < curr_num_edges; ++j) {
 
        size_t curr_edge = curr_cell_to_edge[j];
        size_t * curr_edge_to_cell = edge_to_cell[curr_edge];
        curr_edge_to_cell += *curr_edge_to_cell != SIZE_MAX;
        YAC_ASSERT(
          *curr_edge_to_cell == SIZE_MAX,
          "ERROR(generate_edge_to_cell): "
          "more than two cells point to a single edge")
        *curr_edge_to_cell = i;
      }
    }
  }
 
  return edge_to_cell;
}
 
static void generate_dist_remote_points(
  struct proc_sphere_part_node * proc_sphere_part,
  struct basic_grid_data * grid,
  struct remote_points ** dist_cell_owner,
  struct remote_points ** dist_vertex_owner,
  struct remote_points ** dist_edge_owner, MPI_Comm comm) {
 
  // generate global ids if they are missing
  generate_global_ids(proc_sphere_part, grid, comm);
 
  // determine dist ranks for all cells based on their bounding circles
  int * dist_cell_ranks;
  int * dist_cell_rank_counts;
  size_t * dist_cell_rank_offsets;
  determine_dist_cell_ranks(
    proc_sphere_part, grid, comm,
    &dist_cell_ranks, &dist_cell_rank_counts, &dist_cell_rank_offsets);
 
  int * rank_buffer = NULL;
  size_t rank_buffer_size = 0;
  size_t rank_buffer_array_size = 0;
 
  int * num_ranks_buffer =
    xrealloc(dist_cell_rank_counts,
             (grid->num_cells + grid->num_vertices + grid->num_edges) *
             sizeof(*num_ranks_buffer));
  int * num_cell_ranks = num_ranks_buffer;
  int * num_vertex_ranks = num_ranks_buffer + grid->num_cells;
  int * num_edge_ranks =
    num_ranks_buffer + grid->num_cells + grid->num_vertices;
 
  for (size_t i = 0; i < grid->num_vertices; ++i) {
    if (grid->core_vertex_mask[i]) {
      int num_cells = grid->num_cells_per_vertex[i];
      size_t * cells =
        grid->vertex_to_cell + grid->vertex_to_cell_offsets[i];
 
      size_t max_num_ranks = 0;
      for (int j = 0; j < num_cells; ++j)
        max_num_ranks += (size_t)(num_cell_ranks[cells[j]]);
 
      ENSURE_ARRAY_SIZE(rank_buffer, rank_buffer_array_size,
                        rank_buffer_size + (size_t)max_num_ranks);
 
      int * curr_ranks = rank_buffer + rank_buffer_size;
      int curr_num_ranks = 0;
 
      for (int j = 0; j < num_cells; ++j) {
 
        size_t curr_cell_idx = cells[j];
        int curr_num_cell_ranks = num_cell_ranks[curr_cell_idx];
        int * curr_cell_ranks =
          dist_cell_ranks + dist_cell_rank_offsets[curr_cell_idx];
 
        for (int k = 0; k < curr_num_cell_ranks; ++k)
          insert_rank(curr_ranks, &curr_num_ranks, curr_cell_ranks[k]);
      }
      rank_buffer_size += (size_t)(num_vertex_ranks[i] = curr_num_ranks);
    } else {
      num_vertex_ranks[i] = 0;
    }
  }
 
  size_t_2_pointer edge_to_cell =
    generate_edge_to_cell(
      grid->cell_to_edge, grid->num_vertices_per_cell,
      grid->core_cell_mask, grid->num_cells, grid->num_edges);
 
  for (size_t i = 0; i < grid->num_edges; ++i) {
    if (grid->core_edge_mask[i]) {
 
      size_t max_num_ranks = 0;
      for (int j = 0; j < 2; ++j)
        if (edge_to_cell[i][j] != SIZE_MAX)
        max_num_ranks += (size_t)(num_cell_ranks[edge_to_cell[i][j]]);
 
      ENSURE_ARRAY_SIZE(rank_buffer, rank_buffer_array_size,
                        rank_buffer_size + (size_t)max_num_ranks);
 
      int * curr_ranks = rank_buffer + rank_buffer_size;
      int curr_num_ranks = 0;
 
      for (int j = 0; j < 2; ++j) {
 
        size_t curr_cell_idx = edge_to_cell[i][j];
        if (curr_cell_idx == SIZE_MAX) continue;
 
        int curr_num_cell_ranks = num_cell_ranks[curr_cell_idx];
        int * curr_cell_ranks =
          dist_cell_ranks + dist_cell_rank_offsets[curr_cell_idx];
 
        for (int k = 0; k < curr_num_cell_ranks; ++k)
          insert_rank(curr_ranks, &curr_num_ranks, curr_cell_ranks[k]);
      }
      rank_buffer_size += (size_t)(num_edge_ranks[i] = curr_num_ranks);
    } else {
      num_edge_ranks[i] = 0;
    }
  }
  free(edge_to_cell);
 
  int comm_size;
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  int * sendcounts, * recvcounts, * sdispls, * rdispls;
  yac_get_comm_buffers(
    3, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
  int * int_buffer =
    xmalloc(4 * (size_t)comm_size * sizeof(*int_buffer));
  int * total_sendcounts = int_buffer + 0 * comm_size;
  int * total_recvcounts = int_buffer + 1 * comm_size;
  int * total_sdispls =    int_buffer + 2 * comm_size;
  int * total_rdispls =    int_buffer + 3 * comm_size;
 
  for (size_t i = 0; i < grid->num_cells; ++i) {
    if (grid->core_cell_mask[i]) {
      int curr_num_ranks = num_cell_ranks[i];
      int * ranks = dist_cell_ranks + dist_cell_rank_offsets[i];
      for (int j = 0; j < curr_num_ranks; ++j) sendcounts[3 * ranks[j] + 0]++;
    }
  }
  size_t offset = 0;
  for (size_t i = 0; i < grid->num_vertices; ++i) {
    if (grid->core_vertex_mask[i]) {
      int curr_num_ranks = num_vertex_ranks[i];
      int * ranks = rank_buffer + offset;
      offset += (size_t)curr_num_ranks;
      for (int j = 0; j < curr_num_ranks; ++j) sendcounts[3 * ranks[j] + 1]++;
    }
  }
  for (size_t i = 0; i < grid->num_edges; ++i) {
    if (grid->core_edge_mask[i]) {
      int curr_num_ranks = num_edge_ranks[i];
      int * ranks = rank_buffer + offset;
      offset += (size_t)curr_num_ranks;
      for (int j = 0; j < curr_num_ranks; ++j) sendcounts[3 * ranks[j] + 2]++;
    }
  }
 
  yac_generate_alltoallv_args(
    3, sendcounts, recvcounts, sdispls, rdispls, comm);
 
  size_t receive_counts[3] = {0,0,0};
  for (int i = 0, saccu = 0, raccu = 0; i < comm_size; ++i) {
    total_sdispls[i] = saccu;
    total_rdispls[i] = raccu;
    total_sendcounts[i] = 0;
    total_recvcounts[i] = 0;
    for (int j = 0; j < 3; ++j) {
      total_sendcounts[i] += sendcounts[3 * i + j];
      total_recvcounts[i] += recvcounts[3 * i + j];
      receive_counts[j] += (size_t)(recvcounts[3 * i + j]);
    }
    saccu += total_sendcounts[i];
    raccu += total_recvcounts[i];
  }
  size_t local_data_count = (size_t)(total_sendcounts[comm_size - 1]) +
                            (size_t)(total_sdispls[comm_size - 1]);
  size_t recv_count = (size_t)(total_recvcounts[comm_size - 1]) +
                      (size_t)(total_rdispls[comm_size - 1]);
 
  struct id_pos * id_pos_buffer =
    xmalloc((local_data_count + recv_count) * sizeof(*id_pos_buffer));
  struct id_pos * id_pos_send_buffer = id_pos_buffer;
  struct id_pos * id_pos_recv_buffer =
    id_pos_buffer + local_data_count;
 
  // pack send buffer
  for (size_t i = 0; i < grid->num_cells; ++i) {
    if (grid->core_cell_mask[i]) {
      int curr_num_ranks = num_cell_ranks[i];
      int * ranks = dist_cell_ranks + dist_cell_rank_offsets[i];
      yac_int global_id = grid->cell_ids[i];
      for (int j = 0; j < curr_num_ranks; ++j) {
        int pos = sdispls[3 * ranks[j] + 0 + 1]++;
        id_pos_send_buffer[pos].global_id = global_id;
        id_pos_send_buffer[pos].orig_pos = i;
      }
    }
  }
  offset = 0;
  for (size_t i = 0; i < grid->num_vertices; ++i) {
    if (grid->core_vertex_mask[i]) {
      int curr_num_ranks = num_vertex_ranks[i];
      int * ranks = rank_buffer + offset;
      offset += (size_t)curr_num_ranks;
      yac_int global_id = grid->vertex_ids[i];
      for (int j = 0; j < curr_num_ranks; ++j) {
        int pos = sdispls[3 * ranks[j] + 1 + 1]++;
        id_pos_send_buffer[pos].global_id = global_id;
        id_pos_send_buffer[pos].orig_pos = i;
      }
    }
  }
  for (size_t i = 0; i < grid->num_edges; ++i) {
    if (grid->core_edge_mask[i]) {
      int curr_num_ranks = num_edge_ranks[i];
      int * ranks = rank_buffer + offset;
      offset += (size_t)curr_num_ranks;
      yac_int global_id = grid->edge_ids[i];
      for (int j = 0; j < curr_num_ranks; ++j) {
        int pos = sdispls[3 * ranks[j] + 2 + 1]++;
        id_pos_send_buffer[pos].global_id = global_id;
        id_pos_send_buffer[pos].orig_pos = i;
      }
    }
  }
  free(num_ranks_buffer);
  free(rank_buffer);
  free(dist_cell_ranks);
  free(dist_cell_rank_offsets);
 
  MPI_Datatype id_pos_dt = yac_get_id_pos_mpi_datatype(comm);
 
  // exchange data
  yac_alltoallv_p2p(
    id_pos_send_buffer, total_sendcounts, total_sdispls,
    id_pos_recv_buffer, total_recvcounts, total_rdispls,
    sizeof(*id_pos_send_buffer), id_pos_dt, comm);
 
  yac_mpi_call(MPI_Type_free(&id_pos_dt), comm);
 
  size_t dist_owner_counts[3] = {0, 0, 0};
  for (int i = 0; i < comm_size; ++i)
    for (int j = 0; j < 3; ++j)
      dist_owner_counts[j] += (size_t)(recvcounts[3 * i + j]);
  size_t max_dist_owner_count =
    MAX(MAX(dist_owner_counts[0], dist_owner_counts[1]), dist_owner_counts[2]);
  struct single_remote_point * temp_buffer =
    xmalloc(max_dist_owner_count * sizeof(*temp_buffer));
 
  struct remote_points * dist_owners[3];
 
  // unpack data
  for (int i = 0; i < 3; ++i) {
 
    size_t count = 0;
    for (int j = 0; j < comm_size; ++j) {
      int curr_recvcount = recvcounts[3 * j + i];
      struct id_pos * curr_id_pos =
        id_pos_recv_buffer + rdispls[3 * j + i];
      for (int k = 0; k < curr_recvcount; ++k, ++count) {
        temp_buffer[count].global_id = curr_id_pos[k].global_id;
        temp_buffer[count].data.orig_pos = curr_id_pos[k].orig_pos;
        temp_buffer[count].data.rank = j;
      }
    }
 
    // sort received global ids
    qsort(temp_buffer, count, sizeof(*temp_buffer),
          compare_single_remote_point_global_id);
 
    struct remote_point * unique_ids = xmalloc(count * sizeof(*unique_ids));
    size_t num_unique_ids = 0;
 
    // determine unique global ids
    yac_int prev_id = (count > 0)?temp_buffer[0].global_id - 1:-1;
    for (size_t j = 0; j < count; ++j) {
 
      yac_int curr_id = temp_buffer[j].global_id;
      if (curr_id != prev_id) {
        prev_id = curr_id;
        unique_ids[num_unique_ids].global_id = curr_id;
        unique_ids[num_unique_ids].data.count = 1;
        num_unique_ids++;
      } else {
        unique_ids[num_unique_ids-1].data.count++;
      }
    }
 
    size_t remote_point_info_buffer_size = 0;
    for (size_t j = 0; j < num_unique_ids; ++j)
      if (unique_ids[j].data.count > 1)
        remote_point_info_buffer_size +=
          (size_t)(unique_ids[j].data.count);
 
    dist_owners[i] =
      xmalloc(remote_point_info_buffer_size * sizeof(struct remote_point_info) +
              sizeof(*(dist_owners[0])));
    dist_owners[i]->data =
      (unique_ids = xrealloc(unique_ids, num_unique_ids * sizeof(*unique_ids)));
    dist_owners[i]->count = num_unique_ids;
 
    for (size_t j = 0, l = 0, offset = 0; j < num_unique_ids; ++j) {
      int curr_count = unique_ids[j].data.count;
      if (curr_count == 1) {
        unique_ids[j].data.data.single = temp_buffer[l].data;
        ++l;
      } else {
        unique_ids[j].data.data.multi = &(dist_owners[i]->buffer[offset]);
        for (int k = 0; k < curr_count; ++k, ++l, ++offset) {
          dist_owners[i]->buffer[offset] = temp_buffer[l].data;
        }
      }
    }
 
    qsort(
      unique_ids, num_unique_ids, sizeof(*unique_ids), compare_remote_point);
  }
 
  free(id_pos_buffer);
  free(temp_buffer);
  free(int_buffer);
  yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
 
  *dist_cell_owner = dist_owners[0];
  *dist_vertex_owner = dist_owners[1];
  *dist_edge_owner = dist_owners[2];
}
 
static void generate_sorted_ids(
  yac_int * ids, size_t count, yac_int ** sorted_ids, size_t ** reorder_idx) {
 
  yac_int * sorted_ids_ = xrealloc(*sorted_ids, count * sizeof(*sorted_ids_));
  size_t * reorder_idx_ = xrealloc(*reorder_idx, count * sizeof(*reorder_idx_));
 
  memcpy(sorted_ids_, ids, count * sizeof(*sorted_ids_));
  for (size_t i = 0; i < count; ++i) reorder_idx_[i] = i;
 
  yac_quicksort_index_yac_int_size_t(sorted_ids_, count, reorder_idx_);
 
  *sorted_ids = sorted_ids_;
  *reorder_idx = reorder_idx_;
}
 
static Xt_xmap generate_xmap_data(
  void * data, size_t count, MPI_Comm comm,
  void(*set_sendcounts)(void*,size_t,int*),
  void(*pack)(void*,size_t,int*,int*,int*)) {
 
  int comm_size;
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  int * sendcounts, * recvcounts, * sdispls, * rdispls;
  yac_get_comm_buffers(
    1, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
 
  set_sendcounts(data, count, sendcounts);
 
  yac_generate_alltoallv_args(
    1, sendcounts, recvcounts, sdispls, rdispls, comm);
  int num_src_msg = 0, num_dst_msg = 0;
  for (int i = 0; i < comm_size; ++i) {
    num_src_msg += (recvcounts[i] > 0);
    num_dst_msg += (sendcounts[i] > 0);
  }
 
  size_t recv_count =
    (size_t)(rdispls[comm_size-1]) + (size_t)(recvcounts[comm_size-1]);
 
  int * pos_buffer =
    xmalloc((recv_count + 2 * count) * sizeof(*pos_buffer));
  int * src_pos_buffer = pos_buffer;
  int * dst_pos_buffer = pos_buffer + recv_count;
  int * send_pos_buffer = pos_buffer + recv_count + count;
 
  // pack send buffer
  pack(data, count, sdispls, dst_pos_buffer, send_pos_buffer);
 
  // redistribute positions of requested data
  yac_alltoallv_int_p2p(
    send_pos_buffer, sendcounts, sdispls,
    src_pos_buffer, recvcounts, rdispls, comm);
 
  struct Xt_com_pos * com_pos =
    xmalloc(((size_t)num_src_msg + (size_t)num_dst_msg) * sizeof(*com_pos));
  struct Xt_com_pos * src_com = com_pos;
  struct Xt_com_pos * dst_com = com_pos + num_src_msg;
 
  // set transfer_pos pointers and transfer_pos counts in com_pos's
  num_src_msg = 0;
  num_dst_msg = 0;
  for (int i = 0; i < comm_size; ++i) {
    if (recvcounts[i] > 0) {
      src_com[num_src_msg].transfer_pos = src_pos_buffer;
      src_com[num_src_msg].num_transfer_pos = recvcounts[i];
      src_com[num_src_msg].rank = i;
      src_pos_buffer += recvcounts[i];
      ++num_src_msg;
    }
    if (sendcounts[i] > 0) {
      dst_com[num_dst_msg].transfer_pos = dst_pos_buffer;
      dst_com[num_dst_msg].num_transfer_pos = sendcounts[i];
      dst_com[num_dst_msg].rank = i;
      dst_pos_buffer += sendcounts[i];
      ++num_dst_msg;
    }
  }
  yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
 
  Xt_xmap xmap =
    xt_xmap_intersection_pos_new(
      num_src_msg, src_com, num_dst_msg, dst_com, comm);
 
  free(com_pos);
  free(pos_buffer);
 
  return xmap;
}
 
static void generate_xmap_set_sendcounts_remote_points(
  void * data, size_t count, int * sendcounts) {
 
  struct remote_point * remote_point_data = (struct remote_point *)data;
 
  for (size_t i = 0; i < count; ++i) {
    struct remote_point_info * curr_info =
      (remote_point_data[i].data.count > 1)?
        (remote_point_data[i].data.data.multi):
        (&(remote_point_data[i].data.data.single));
    sendcounts[curr_info->rank]++;
  }
}
 
static void generate_xmap_pack_remote_points(
  void * data, size_t count, int * sdispls,
  int * dst_pos_buffer, int * send_pos_buffer) {
 
  struct remote_point * remote_point_data = (struct remote_point *)data;
 
  for (size_t i = 0; i < count; ++i) {
    struct remote_point_info * curr_info =
      (remote_point_data[i].data.count > 1)?
        (remote_point_data[i].data.data.multi):
        (&(remote_point_data[i].data.data.single));
    int pos = sdispls[curr_info->rank+1]++;
    dst_pos_buffer[pos] = (int)i;
    send_pos_buffer[pos] = (int)(curr_info->orig_pos);
  }
}
 
static Xt_xmap generate_xmap_from_remote_points(
  struct remote_point * remote_point_data, size_t count, MPI_Comm comm) {
 
  return
    generate_xmap_data(remote_point_data, count, comm,
                       generate_xmap_set_sendcounts_remote_points,
                       generate_xmap_pack_remote_points);
}
 
static void generate_xmap_set_sendcounts_remote_point_info(
  void * data, size_t count, int * sendcounts) {
 
  struct remote_point_info * point_infos = (struct remote_point_info *)data;
 
  for (size_t i = 0; i < count; ++i) sendcounts[point_infos[i].rank]++;
}
 
static void generate_xmap_pack_remote_point_info(
  void * data, size_t count, int * sdispls,
  int * dst_pos_buffer, int * send_pos_buffer) {
 
  struct remote_point_info * point_infos = (struct remote_point_info *)data;
 
  for (size_t i = 0; i < count; ++i) {
    int pos = sdispls[point_infos[i].rank+1]++;
    dst_pos_buffer[pos] = (int)i;
    send_pos_buffer[pos] = (int)(point_infos[i].orig_pos);
  }
}
 
static Xt_xmap generate_xmap_from_remote_point_info(
  struct remote_point_info * point_infos, size_t count, MPI_Comm comm) {
 
  return
    generate_xmap_data(point_infos, count, comm,
                       generate_xmap_set_sendcounts_remote_point_info,
                       generate_xmap_pack_remote_point_info);
}
 
static int compare_single_remote_point_reorder_global_id(
  const void * a, const void * b) {
 
  return (((const struct single_remote_point_reorder *)a)->data.global_id >
          ((const struct single_remote_point_reorder *)b)->data.global_id) -
         (((const struct single_remote_point_reorder *)a)->data.global_id <
          ((const struct single_remote_point_reorder *)b)->data.global_id);
}
 
static int compare_single_remote_point_reorder_reorder_idx(
  const void * a, const void * b) {
 
  return (((const struct single_remote_point_reorder *)a)->reorder_idx >
          ((const struct single_remote_point_reorder *)b)->reorder_idx) -
         (((const struct single_remote_point_reorder *)a)->reorder_idx <
          ((const struct single_remote_point_reorder *)b)->reorder_idx);
}
 
static void compute_basic_ve_data(
  size_t max_num_cell_ve, size_t num_cells, int * num_ve_per_cell,
  struct single_remote_point_reorder * cell_ve_point_data,
  yac_int ** ids, size_t * num_ids, size_t ** cell_to_ve, Xt_xmap * xmap,
  MPI_Comm comm) {
 
  size_t num_total_ve = num_cells * max_num_cell_ve;
 
  size_t total_num_cell_ve = 0;
  for (size_t i = 0, k = 0, l = 0; i < num_cells; ++i) {
    size_t num_vertices = (size_t)(num_ve_per_cell[i]);
    total_num_cell_ve += num_vertices;
    size_t j = 0;
    for (; j < num_vertices; ++j, ++k, ++l) {
      cell_ve_point_data[k].reorder_idx = l;
    }
    for (; j < max_num_cell_ve; ++j, ++k) {
      cell_ve_point_data[k].data.global_id = XT_INT_MAX;
      cell_ve_point_data[k].reorder_idx = SIZE_MAX;
    }
  }
  qsort(cell_ve_point_data, num_total_ve, sizeof(*cell_ve_point_data),
        compare_single_remote_point_reorder_global_id);
  yac_int * ids_ = xmalloc(total_num_cell_ve * sizeof(*ids_));
  size_t num_ids_ = 0;
  size_t * cell_to_ve_ =
    xmalloc(total_num_cell_ve * sizeof(*cell_to_ve_));
  struct remote_point_info * ve_point_info =
    xmalloc(total_num_cell_ve * sizeof(*ve_point_info));
  yac_int prev_global_id = XT_INT_MAX;
  for (size_t i = 0; i < num_total_ve; ++i) {
    yac_int curr_global_id = cell_ve_point_data[i].data.global_id;
    size_t curr_reorder_idx = cell_ve_point_data[i].reorder_idx;
    if (curr_global_id == XT_INT_MAX) break;
    if (curr_global_id != prev_global_id) {
      ids_[num_ids_] = (prev_global_id = curr_global_id);
      ve_point_info[num_ids_] = cell_ve_point_data[i].data.data;
      num_ids_++;
    }
    cell_to_ve_[curr_reorder_idx] = num_ids_ - 1;
  }
 
  *ids = xrealloc(ids_, num_ids_ * sizeof(*ids_));
  *num_ids = num_ids_;
  *cell_to_ve = cell_to_ve_;
  *xmap = generate_xmap_from_remote_point_info(ve_point_info, num_ids_, comm);
  free(ve_point_info);
}
 
static MPI_Datatype yac_get_coordinate_mpi_datatype(MPI_Comm comm) {
 
  MPI_Datatype coord_dt;
  yac_mpi_call(MPI_Type_contiguous(3, MPI_DOUBLE, &coord_dt), comm);
  yac_mpi_call(MPI_Type_commit(&coord_dt), comm);
  return coord_dt;
}
 
static struct yac_field_data field_data_init(
  struct yac_field_data * orig_field_data, size_t dist_size,
  Xt_redist redist_mask, Xt_redist redist_coords, MPI_Comm comm) {
 
  struct yac_field_data dist_field_data = yac_field_data_init();
 
  uint64_t counts[2], max_counts[2];
  if (orig_field_data != NULL) {
    counts[0] = orig_field_data->masks_count;
    counts[1] = orig_field_data->coordinates_count;
  } else {
    counts[0] = 0;
    counts[1] = 0;
  }
  yac_mpi_call(
    MPI_Allreduce(
      counts, max_counts, 2, MPI_UINT64_T, MPI_MAX, comm), comm);
  YAC_ASSERT(
    (orig_field_data == NULL) ||
    ((counts[0] == max_counts[0]) && (counts[1] == max_counts[1])),
    "ERROR(field_data_init): inconsistent number of masks or coordinates")
 
  int * data_available_flag =
    xcalloc(max_counts[0] + max_counts[1], sizeof(*data_available_flag));
 
  if (orig_field_data != NULL) {
    for (size_t i = 0; i < orig_field_data->masks_count; ++i)
      data_available_flag[i] = orig_field_data->masks[i] != NULL;
    for (size_t i = 0; i < orig_field_data->coordinates_count; ++i)
      data_available_flag[i + orig_field_data->masks_count] =
        orig_field_data->coordinates[i] != NULL;
  }
 
  yac_mpi_call(
    MPI_Allreduce(
      MPI_IN_PLACE, data_available_flag,
      (int)(max_counts[0] + max_counts[1]), MPI_INT, MPI_MAX, comm), comm);
 
  if (orig_field_data != NULL) {
    for (size_t i = 0; i < orig_field_data->masks_count; ++i)
      YAC_ASSERT(
        data_available_flag[i] == (orig_field_data->masks[i] != NULL),
        "ERROR(field_data_init): inconsistent availability of masks")
 
    for (size_t i = 0; i < orig_field_data->coordinates_count; ++i)
      YAC_ASSERT(
        data_available_flag[i + orig_field_data->masks_count] ==
        (orig_field_data->coordinates[i] != NULL),
        "ERROR(field_data_init): inconsistent availability of coordinates")
  }
 
  for (uint64_t i = 0; i < max_counts[0]; ++i) {
    int * dist_mask = NULL;
    if (data_available_flag[i]) {
      dist_mask = xmalloc(dist_size * sizeof(*dist_mask));
      int * orig_mask = (orig_field_data != NULL)?orig_field_data->masks[i]:NULL;
      xt_redist_s_exchange1(redist_mask, orig_mask, dist_mask);
    }
    yac_field_data_add_mask_nocpy(&dist_field_data, dist_mask);
  }
 
  for (uint64_t i = 0; i < max_counts[1]; ++i) {
    coordinate_pointer dist_coordinates = NULL;
    if (data_available_flag[i + max_counts[0]]) {
      dist_coordinates = xmalloc(dist_size * sizeof(*dist_coordinates));
      coordinate_pointer orig_coordinates =
        (orig_field_data != NULL)?orig_field_data->coordinates[i]:NULL;
      xt_redist_s_exchange1(redist_coords, orig_coordinates, dist_coordinates);
    }
    yac_field_data_add_coordinates_nocpy(&dist_field_data, dist_coordinates);
  }
 
  free(data_available_flag);
 
  return dist_field_data;
}
 
static void generate_dist_grid(
  struct dist_grid * dist_grid, struct proc_sphere_part_node * proc_sphere_part,
  struct yac_basic_grid * grid, MPI_Comm comm) {
 
  struct basic_grid_data * grid_data = yac_basic_grid_get_data(grid);
 
  int comm_rank;
  yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
 
  int max_num_vertices_per_cell = 0;
  { // compute the global maximum number of vertices per cell
    for (size_t i = 0; i < grid_data->num_cells; ++i)
      if (grid_data->num_vertices_per_cell[i] > max_num_vertices_per_cell)
        max_num_vertices_per_cell = grid_data->num_vertices_per_cell[i];
    yac_mpi_call(MPI_Allreduce(MPI_IN_PLACE, &max_num_vertices_per_cell, 1,
                               MPI_INT, MPI_MAX, comm), comm);
  }
 
  struct remote_points * dist_cell_owner,
                       * dist_vertex_owner,
                       * dist_edge_owner;
 
  generate_dist_remote_points(
    proc_sphere_part, grid_data,
    &dist_cell_owner, &dist_vertex_owner, &dist_edge_owner, comm);
 
  Xt_xmap xmap_cell_orig_to_dist =
    generate_xmap_from_remote_points(
      dist_cell_owner->data, dist_cell_owner->count, comm);
 
  MPI_Datatype dt_single_remote_point =
    yac_get_n_single_remote_point_reorder_mpi_datatype(
      max_num_vertices_per_cell, comm);
  MPI_Datatype dt_coord = yac_get_coordinate_mpi_datatype(comm);
 
  Xt_redist redist_cell_yac_int,
            redist_cell_int,
            redist_cell_corner_point_data,
            redist_cell_coords;
  redist_cell_yac_int = xt_redist_p2p_new(xmap_cell_orig_to_dist, Xt_int_dt);
  redist_cell_int = xt_redist_p2p_new(xmap_cell_orig_to_dist, MPI_INT);
  redist_cell_corner_point_data =
    xt_redist_p2p_new(xmap_cell_orig_to_dist, dt_single_remote_point);
  yac_mpi_call(MPI_Type_free(&dt_single_remote_point), comm);
  redist_cell_coords = xt_redist_p2p_new(xmap_cell_orig_to_dist, dt_coord);
 
  xt_xmap_delete(xmap_cell_orig_to_dist);
 
  size_t num_cells = dist_cell_owner->count;
  size_t max_num_cell_corners =
    (size_t)max_num_vertices_per_cell * num_cells;
  size_t max_num_local_cell_corners =
    (size_t)max_num_vertices_per_cell * grid_data->num_cells;
  yac_int * cell_ids = xmalloc(num_cells * sizeof(*cell_ids));
  int * num_vertices_per_cell =
    xmalloc(num_cells * sizeof(*num_vertices_per_cell));
  struct single_remote_point_reorder *
    cell_corner_point_data =
    xmalloc(2 * (max_num_cell_corners + max_num_local_cell_corners) *
            sizeof(*cell_corner_point_data));
  struct single_remote_point_reorder * cell_edge_point_data =
    cell_corner_point_data + max_num_cell_corners;
  struct single_remote_point_reorder * local_cell_corner_point_data =
    cell_corner_point_data + 2 * max_num_cell_corners;
  struct single_remote_point_reorder * local_cell_edge_point_data =
    cell_corner_point_data + 2 * max_num_cell_corners +
    max_num_local_cell_corners;
  for (size_t i = 0; i < grid_data->num_cells; ++i){
    int num_corners = (size_t)(grid_data->num_vertices_per_cell[i]);
    size_t * corners = grid_data->cell_to_vertex + grid_data->cell_to_vertex_offsets[i];
    size_t * edges = grid_data->cell_to_edge + grid_data->cell_to_edge_offsets[i];
    struct single_remote_point_reorder * curr_local_cell_corner_point_data =
      local_cell_corner_point_data + i * (size_t)max_num_vertices_per_cell;
    struct single_remote_point_reorder * curr_local_cell_edge_point_data =
      local_cell_edge_point_data + i * (size_t)max_num_vertices_per_cell;
    for (size_t j = 0; j < num_corners; ++j) {
      curr_local_cell_corner_point_data[j].data.global_id = grid_data->vertex_ids[corners[j]];
      curr_local_cell_corner_point_data[j].data.data.rank = comm_rank;
      curr_local_cell_corner_point_data[j].data.data.orig_pos = corners[j];
      curr_local_cell_edge_point_data[j].data.global_id = grid_data->edge_ids[edges[j]];
      curr_local_cell_edge_point_data[j].data.data.rank = comm_rank;
      curr_local_cell_edge_point_data[j].data.data.orig_pos = edges[j];
    }
  }
 
  xt_redist_s_exchange1(redist_cell_yac_int, grid_data->cell_ids, cell_ids);
  xt_redist_s_exchange1(
    redist_cell_int, grid_data->num_vertices_per_cell, num_vertices_per_cell);
  xt_redist_s_exchange1(
    redist_cell_corner_point_data,
    local_cell_corner_point_data, cell_corner_point_data);
  xt_redist_s_exchange1(
    redist_cell_corner_point_data,
    local_cell_edge_point_data, cell_edge_point_data);
 
  struct yac_field_data cell_field_data =
    field_data_init(
      yac_basic_grid_get_field_data(grid, CELL), num_cells,
      redist_cell_int, redist_cell_coords, comm);
 
  xt_redist_delete(redist_cell_coords);
  xt_redist_delete(redist_cell_corner_point_data);
  xt_redist_delete(redist_cell_yac_int);
  xt_redist_delete(redist_cell_int);
 
  yac_int * vertex_ids;
  size_t num_vertices;
  size_t * cell_to_vertex;
  size_t * cell_to_vertex_offsets =
    xmalloc(num_cells * sizeof(*cell_to_vertex_offsets));
  Xt_xmap xmap_vertex_orig_to_dist;
  for (size_t i = 0, accu = 0; i < num_cells; ++i) {
    cell_to_vertex_offsets[i] = accu;
    accu += (size_t)(num_vertices_per_cell[i]);
  }
  compute_basic_ve_data(
    max_num_vertices_per_cell, num_cells, num_vertices_per_cell,
    cell_corner_point_data, &vertex_ids, &num_vertices,
    &cell_to_vertex, &xmap_vertex_orig_to_dist, comm);
 
  yac_int * edge_ids;
  size_t num_edges;
  size_t * cell_to_edge;
  size_t * cell_to_edge_offsets = cell_to_vertex_offsets;
  Xt_xmap xmap_edge_orig_to_dist;
  compute_basic_ve_data(
    max_num_vertices_per_cell, num_cells, num_vertices_per_cell,
    cell_edge_point_data, &edge_ids, &num_edges,
    &cell_to_edge, &xmap_edge_orig_to_dist, comm);
 
  free(cell_corner_point_data);
 
  // generate redists for redistribution of vertex data
  Xt_redist redist_vertex_coords =
    xt_redist_p2p_new(xmap_vertex_orig_to_dist, dt_coord);
  Xt_redist redist_vertex_int =
    xt_redist_p2p_new(xmap_vertex_orig_to_dist, MPI_INT);
  xt_xmap_delete(xmap_vertex_orig_to_dist);
 
  // get vertex coordinates
  coordinate_pointer vertex_coordinates =
    xmalloc(num_vertices * sizeof(*vertex_coordinates));
  xt_redist_s_exchange1(
    redist_vertex_coords, grid_data->vertex_coordinates, vertex_coordinates);
 
  struct yac_field_data vertex_field_data =
    field_data_init(
      yac_basic_grid_get_field_data(grid, CORNER), num_vertices,
      redist_vertex_int, redist_vertex_coords, comm);
 
  xt_redist_delete(redist_vertex_int);
  xt_redist_delete(redist_vertex_coords);
 
  Xt_redist redist_edge_int,
            redist_edge_coords,
            redist_edge_2yac_int;
  { // generate redists for redistribution of edge data
 
    MPI_Datatype dt_2yac_int;
    yac_mpi_call(MPI_Type_contiguous(2, Xt_int_dt, &dt_2yac_int), comm);
    yac_mpi_call(MPI_Type_commit(&dt_2yac_int), comm);
 
    redist_edge_int = xt_redist_p2p_new(xmap_edge_orig_to_dist, MPI_INT);
    redist_edge_coords = xt_redist_p2p_new(xmap_edge_orig_to_dist, dt_coord);
    redist_edge_2yac_int = xt_redist_p2p_new(xmap_edge_orig_to_dist, dt_2yac_int);
 
    yac_mpi_call(MPI_Type_free(&dt_2yac_int), comm);
 
  }
  xt_xmap_delete(xmap_edge_orig_to_dist);
 
  enum yac_edge_type * edge_type = xmalloc(num_edges * sizeof(*edge_type));
 
  { // get edge types
    int * temp_edge_type_src, * temp_edge_type_dst;
    if (sizeof(*edge_type) == sizeof(int)) {
      temp_edge_type_src = (int*)(grid_data->edge_type);
      temp_edge_type_dst = (int*)edge_type;
    } else {
        temp_edge_type_src =
          xmalloc(grid_data->num_edges * sizeof(*temp_edge_type_src));
        for (size_t i = 0; i < grid_data->num_edges; ++i)
          temp_edge_type_src[i] = (int)(grid_data->edge_type[i]);
        temp_edge_type_dst = xmalloc(num_edges * sizeof(*temp_edge_type_dst));
        for (size_t i = 0; i < num_edges; ++i)
          temp_edge_type_dst[i] = (int)(edge_type[i]);
    }
 
    xt_redist_s_exchange1(
      redist_edge_int, temp_edge_type_src, temp_edge_type_dst);
 
    if (sizeof(*edge_type) != sizeof(int)) {
 
      for (size_t i = 0; i < num_edges; ++i)
        edge_type[i] = (enum yac_edge_type)(temp_edge_type_dst[i]);
 
      free(temp_edge_type_src);
      free(temp_edge_type_dst);
    }
  }
 
  struct yac_field_data edge_field_data =
    field_data_init(
      yac_basic_grid_get_field_data(grid, EDGE), num_edges,
      redist_edge_int, redist_edge_coords, comm);
 
  xt_redist_delete(redist_edge_int);
  xt_redist_delete(redist_edge_coords);
 
  struct bounding_circle * cell_bnd_circles =
    xmalloc(num_cells * sizeof(*cell_bnd_circles));
  {
    struct grid_cell cell;
    cell.coordinates_xyz = xmalloc((size_t)max_num_vertices_per_cell *
                                   sizeof(*(cell.coordinates_xyz)));
    cell.edge_type = xmalloc((size_t)max_num_vertices_per_cell *
                             sizeof(*(cell.edge_type)));
 
    for (size_t i = 0; i < num_cells; ++i) {
      size_t * curr_cell_to_vertex =
        cell_to_vertex + cell_to_vertex_offsets[i];
      size_t * curr_cell_to_edge =
        cell_to_edge + cell_to_edge_offsets[i];
      for (int j = 0; j < num_vertices_per_cell[i]; ++j) {
        coordinate_pointer curr_vertex_coords =
          vertex_coordinates + curr_cell_to_vertex[j];
        for (int k = 0; k < 3; ++k)
          cell.coordinates_xyz[j][k] = (*curr_vertex_coords)[k];
        cell.edge_type[j] = edge_type[curr_cell_to_edge[j]];
      }
      cell.num_corners = num_vertices_per_cell[i];
      if (cell.num_corners > 0)
        yac_get_cell_bounding_circle(cell, cell_bnd_circles + i);
      else
        cell_bnd_circles[i] =
          (struct bounding_circle) {
            .base_vector = {1.0, 0.0, 0.0},
            .inc_angle = SIN_COS_ZERO,
            .sq_crd = DBL_MAX};
    }
    free(cell.edge_type);
    free(cell.coordinates_xyz);
  }
 
  size_t_2_pointer edge_to_vertex =
    xmalloc(num_edges * sizeof(*edge_to_vertex));
 
  // get edge to vertex
  {
    yac_int * vertex_id_buffer =
      xmalloc(
        2 * (grid_data->num_edges + num_edges) * sizeof(*vertex_id_buffer));
    yac_int * grid_edge_vertex_ids = vertex_id_buffer;
    yac_int * edge_vertex_ids = vertex_id_buffer + 2 * grid_data->num_edges;
 
    size_t * grid_edge_to_vertex = &(grid_data->edge_to_vertex[0][0]);
 
    for (size_t i = 0; i < 2 * grid_data->num_edges; ++i)
      grid_edge_vertex_ids[i] =
        grid_data->vertex_ids[grid_edge_to_vertex[i]];
 
    xt_redist_s_exchange1(
      redist_edge_2yac_int, grid_edge_vertex_ids, edge_vertex_ids);
 
    id2idx(edge_vertex_ids, &(edge_to_vertex[0][0]), 2 * num_edges,
           vertex_ids, num_vertices);
 
    free(vertex_id_buffer);
  }
 
  xt_redist_delete(redist_edge_2yac_int);
 
  size_t num_total_cells = num_cells;
  size_t num_total_vertices = num_vertices;
  size_t num_total_edges = num_edges;
 
  yac_mpi_call(MPI_Type_free(&dt_coord), comm);
 
  dist_grid->comm = comm;
  dist_grid->vertex_coordinates = vertex_coordinates;
  dist_grid->cell_ids = cell_ids;
  dist_grid->vertex_ids = vertex_ids;
  dist_grid->edge_ids = edge_ids;
  dist_grid->num_cells = num_cells;
  dist_grid->num_vertices = num_vertices;
  dist_grid->num_edges = num_edges;
  dist_grid->num_vertices_per_cell = num_vertices_per_cell;
  dist_grid->cell_to_vertex = cell_to_vertex;
  dist_grid->cell_to_vertex_offsets = cell_to_vertex_offsets;
  dist_grid->cell_to_edge = cell_to_edge;
  dist_grid->cell_to_edge_offsets = cell_to_edge_offsets;
  dist_grid->edge_to_vertex = edge_to_vertex;
  dist_grid->cell_bnd_circles = cell_bnd_circles;
  dist_grid->edge_type = edge_type;
  dist_grid->sorted_cell_ids = NULL;
  dist_grid->sorted_vertex_ids = NULL;
  dist_grid->sorted_edge_ids = NULL;
  dist_grid->sorted_cell_reorder_idx = NULL;
  dist_grid->sorted_vertex_reorder_idx = NULL;
  dist_grid->sorted_edge_reorder_idx = NULL;
  generate_sorted_ids(cell_ids, num_cells, &(dist_grid->sorted_cell_ids),
                      &(dist_grid->sorted_cell_reorder_idx));
  generate_sorted_ids(vertex_ids, num_vertices, &(dist_grid->sorted_vertex_ids),
                      &(dist_grid->sorted_vertex_reorder_idx));
  generate_sorted_ids(edge_ids, num_edges, &(dist_grid->sorted_edge_ids),
                      &(dist_grid->sorted_edge_reorder_idx));
  dist_grid->core_vertex_mask =
    determine_core_vertex_mask(
      proc_sphere_part, comm_rank, vertex_coordinates, num_total_vertices);
  dist_grid->core_cell_mask =
    determine_cell_core_mask(
      dist_grid->core_vertex_mask, dist_grid, comm_rank == 0);
  dist_grid->core_edge_mask =
    determine_edge_core_mask(dist_grid->core_vertex_mask, dist_grid);
  dist_grid->num_total_cells = num_total_cells;
  dist_grid->num_total_vertices = num_total_vertices;
  dist_grid->num_total_edges = num_total_edges;
  dist_grid->field_data.cell = cell_field_data;
  dist_grid->field_data.vertex = vertex_field_data;
  dist_grid->field_data.edge = edge_field_data;
 
  generate_dist_grid_remote_point_infos(
    dist_grid, proc_sphere_part,
    dist_cell_owner, dist_vertex_owner, dist_edge_owner);
 
  free(dist_cell_owner->data);
  free(dist_cell_owner);
  free(dist_vertex_owner->data);
  free(dist_vertex_owner);
  free(dist_edge_owner->data);
  free(dist_edge_owner);
}
 
static void setup_search_data(struct dist_grid_pair * dist_grid_pair) {
 
  // build sphere part for vertices
  for (int i = 0; i < 2; ++i)
    dist_grid_pair->vertex_sphere_part[i] =
      yac_point_sphere_part_search_new(
        dist_grid_pair->dist_grid[i].num_vertices,
        dist_grid_pair->dist_grid[i].vertex_coordinates,
        dist_grid_pair->dist_grid[i].vertex_ids);
 
  // build sphere part for bounding circle of cells
  for (int i = 0; i < 2; ++i) {
    dist_grid_pair->cell_sphere_part[i] =
      yac_bnd_sphere_part_search_new(
        dist_grid_pair->dist_grid[i].cell_bnd_circles,
        dist_grid_pair->dist_grid[i].num_cells);
  }
}
 
static void get_grid_cells(
  struct yac_basic_grid * grid, struct dist_cell * cells, MPI_Comm comm) {
 
  struct basic_grid_data * grid_data = yac_basic_grid_get_data(grid);
 
  size_t num_cells = grid_data->num_cells;
  for (size_t i = 0; i < num_cells; ++i) {
 
    double * cell_coord = cells[i].coord;
    for (int j = 0; j < 3; ++j) cell_coord[j] = 0.0;
 
    if (grid_data->num_vertices_per_cell[i] == 0) continue;
 
    size_t * curr_cell_to_vertex =
      grid_data->cell_to_vertex + grid_data->cell_to_vertex_offsets[i];
    // compute the middle point of the cell
    for (int j = 0; j < grid_data->num_vertices_per_cell[i]; ++j) {
      coordinate_pointer curr_vertex_coords =
        grid_data->vertex_coordinates + curr_cell_to_vertex[j];
      for (int k = 0; k < 3; ++k)
        cell_coord[k] += (*curr_vertex_coords)[k];
    }
    normalise_vector(cell_coord);
  }
}
 
struct dist_grid_pair * yac_dist_grid_pair_new(
  struct yac_basic_grid * grid_a, int grid_id_a,
  struct yac_basic_grid * grid_b, int grid_id_b, MPI_Comm comm) {
 
  YAC_ASSERT(
    grid_id_a != grid_id_b,"ERROR(yac_dist_grid_pair_new): identical grid ids")
 
  MPI_Comm comm_copy;
  yac_mpi_call(MPI_Comm_dup(comm, &comm_copy), comm);
  comm = comm_copy;
 
  if (grid_id_a > grid_id_b) {
    int grid_id_swap = grid_id_a;
    grid_id_a = grid_id_b;
    grid_id_b = grid_id_swap;
    struct yac_basic_grid * grid_swap = grid_a;
    grid_a = grid_b;
    grid_b = grid_swap;
  }
 
  struct dist_grid_pair * grid_pair = xmalloc(1 * sizeof(*grid_pair));
 
  size_t num_cells_a = yac_basic_grid_get_data(grid_a)->num_cells;
  size_t num_cells_b = yac_basic_grid_get_data(grid_b)->num_cells;
 
  size_t num_cells = num_cells_a + num_cells_b;
  struct dist_cell * cells = xmalloc(num_cells * sizeof(*cells));
  get_grid_cells(grid_a, cells, comm);
  get_grid_cells(grid_b, cells + num_cells_a, comm);
 
  grid_pair->grid_ids[0] = grid_id_a;
  grid_pair->grid_ids[1] = grid_id_b;
  grid_pair->comm = comm;
 
  // redistribute all points and build parallel sphere part for vertices
  grid_pair->proc_sphere_part =
    yac_redistribute_cells(&cells, &num_cells, comm);
  free(cells);
 
  // generate dist_grids
  generate_dist_grid(
    &(grid_pair->dist_grid[0]), grid_pair->proc_sphere_part, grid_a, comm);
  generate_dist_grid(
    &(grid_pair->dist_grid[1]), grid_pair->proc_sphere_part, grid_b, comm);
 
  // build search data structures for local cells and vertices in
  // distributed grid
  setup_search_data(grid_pair);
 
  return grid_pair;
}
 
MPI_Comm yac_dist_grid_pair_get_MPI_Comm(struct dist_grid_pair * grid_pair) {
  return grid_pair->comm;
}
 
struct dist_grid * yac_dist_grid_pair_get_dist_grid(
  struct dist_grid_pair * grid_pair, int grid_id) {
 
  struct dist_grid * grid = NULL;
  for (int i = 0; (i < 2) && (grid == NULL); ++i)
    if (grid_id == grid_pair->grid_ids[i])
      grid = &(grid_pair->dist_grid[i]);
  YAC_ASSERT(
    grid != NULL, "ERROR(yac_dist_grid_pair_get_dist_grid): invalid grid_id")
  return grid;
}
 
struct const_basic_grid_data * yac_dist_grid_get_basic_grid_data(
  struct dist_grid * dist_grid) {
 
  return (struct const_basic_grid_data *)dist_grid;
}
 
size_t yac_dist_grid_get_local_count(
  struct dist_grid * dist_grid, enum yac_location location) {
 
  YAC_ASSERT(
    (location == CELL) || (location == CORNER) || (location == EDGE),
    "ERROR(yac_dist_grid_get_local_count): invalid location")
  switch (location) {
    default:
    case(CELL):
      return
        determine_local_count(
          dist_grid->num_cells, dist_grid->core_cell_mask);
    case (CORNER):
      return
        determine_local_count(
          dist_grid->num_vertices, dist_grid->core_vertex_mask);
    case (EDGE):
      return
        determine_local_count(
          dist_grid->num_edges, dist_grid->core_edge_mask);
  };
}
 
static size_t yac_dist_grid_get_count(
  struct dist_grid * dist_grid, enum yac_location location) {
 
  YAC_ASSERT(
    (location == CELL) || (location == CORNER) || (location == EDGE),
    "ERROR(yac_dist_grid_get_count): invalid location")
  switch (location) {
    default:
    case(CELL): return dist_grid->num_total_cells;
    case (CORNER): return dist_grid->num_total_vertices;
    case (EDGE): return dist_grid->num_total_edges;
  };
}
 
static int const * yac_dist_grid_get_core_mask(
  struct dist_grid * dist_grid, enum yac_location location) {
 
  YAC_ASSERT(
    (location == CELL) || (location == CORNER) || (location == EDGE),
    "ERROR(yac_dist_grid_get_core_mask): invalid location")
  switch (location) {
    default:
    case(CELL): return dist_grid->core_cell_mask;
    case (CORNER): return dist_grid->core_vertex_mask;
    case (EDGE): return dist_grid->core_edge_mask;
  };
}
 
static yac_int const * yac_dist_grid_get_global_ids(
  struct dist_grid * dist_grid, enum yac_location location) {
 
  YAC_ASSERT(
    (location == CELL) || (location == CORNER) || (location == EDGE),
    "ERROR(yac_dist_grid_get_global_ids): invalid location")
  switch (location) {
    default:
    case(CELL): return dist_grid->cell_ids;
    case (CORNER): return dist_grid->vertex_ids;
    case (EDGE): return dist_grid->edge_ids;
  };
}
 
void yac_dist_grid_get_local_unmasked_points(
  struct dist_grid * dist_grid, struct interp_field field,
  size_t ** indices, size_t * count) {
 
  int const * mask = yac_dist_grid_get_field_mask(dist_grid, field);
 
  size_t total_count = yac_dist_grid_get_count(dist_grid, field.location);
  int const * core_mask =
    yac_dist_grid_get_core_mask(dist_grid, field.location);
 
  size_t * temp_indices = xmalloc(total_count * sizeof(*temp_indices));
 
  size_t count_ = 0;
 
  if (mask != NULL) {
    for (size_t i = 0; i < total_count; ++i)
      if (core_mask[i] && mask[i]) temp_indices[count_++] = i;
  } else {
    for (size_t i = 0; i < total_count; ++i)
      if (core_mask[i]) temp_indices[count_++] = i;
  }
 
  *indices = xrealloc(temp_indices, count_ * sizeof(**indices));
  *count = count_;
}
 
static struct yac_field_data * yac_dist_grid_get_field_data(
  struct dist_grid * dist_grid, enum yac_location location) {
 
  YAC_ASSERT(
    (location == CELL) || (location == CORNER) || (location == EDGE),
    "ERROR(yac_dist_grid_get_field_data): invalid location")
 
  switch (location) {
    default:
    case(CELL): return &(dist_grid->field_data.cell);
    case(CORNER): return &(dist_grid->field_data.vertex);
    case(EDGE): return &(dist_grid->field_data.edge);
  };
}
 
int const * yac_dist_grid_get_field_mask(
  struct dist_grid * dist_grid, struct interp_field field) {
 
  if (field.masks_idx == SIZE_MAX) return NULL;
 
  struct yac_field_data * data =
    yac_dist_grid_get_field_data(dist_grid, field.location);
 
  YAC_ASSERT(
    field.masks_idx < data->masks_count,
    "ERROR(yac_dist_grid_get_field_mask): invalid mask_idx")
 
  return data->masks[field.masks_idx];
}
 
coordinate_pointer yac_dist_grid_get_field_coords(
  struct dist_grid * dist_grid, struct interp_field field) {
 
  coordinate_pointer coords = NULL;
 
  if (field.coordinates_idx != SIZE_MAX) {
 
    struct yac_field_data * data =
      yac_dist_grid_get_field_data(dist_grid, field.location);
 
    YAC_ASSERT(
      field.coordinates_idx < data->coordinates_count,
      "ERROR(yac_dist_grid_get_field_coords): invalid coord_idx")
 
    coords = data->coordinates[field.coordinates_idx];
  }
 
  // if no field coordinates are defined, but the location is at the corners of
  // of the grid cells, return coordinates of them
  if ((coords == NULL) && (field.location == CORNER))
    coords = dist_grid->vertex_coordinates;
 
  return coords;
}
 
size_t yac_dist_grid_get_unmasked_local_count(
  struct dist_grid * dist_grid, struct interp_field field) {
 
  size_t total_count = yac_dist_grid_get_count(dist_grid, field.location);
 
  int const * mask = yac_dist_grid_get_field_mask(dist_grid, field);
  if (mask == NULL)
    return yac_dist_grid_get_local_count(dist_grid, field.location);
 
  int const * core_mask =
    yac_dist_grid_get_core_mask(dist_grid, field.location);
 
  size_t unmasked_local_count = 0;
  for (size_t i = 0; i < total_count; ++i)
    if (core_mask[i] && mask[i]) ++unmasked_local_count;
  return unmasked_local_count;
}
 
static void yac_remote_point_infos_free(
  struct remote_point_infos * point_infos, size_t count) {
 
  for (size_t i = 0; i < count; ++i)
    if (point_infos[i].count > 1) free(point_infos[i].data.multi);
  free(point_infos);
}
 
static void yac_dist_grid_free(struct dist_grid grid) {
 
  free(grid.vertex_coordinates);
  free(grid.cell_ids);
  free(grid.vertex_ids);
  free(grid.edge_ids);
  free(grid.num_vertices_per_cell);
  free(grid.cell_to_vertex);
  free(grid.cell_to_vertex_offsets);
  free(grid.cell_to_edge);
  free(grid.cell_bnd_circles);
  free(grid.edge_type);
  free(grid.edge_to_vertex);
  free(grid.core_cell_mask);
  free(grid.core_vertex_mask);
  free(grid.core_edge_mask);
  free(grid.sorted_cell_ids);
  free(grid.sorted_vertex_ids);
  free(grid.sorted_edge_ids);
  free(grid.sorted_cell_reorder_idx);
  free(grid.sorted_vertex_reorder_idx);
  free(grid.sorted_edge_reorder_idx);
  yac_remote_point_infos_free(grid.cell_owners, grid.num_total_cells);
  yac_remote_point_infos_free(grid.vertex_owners, grid.num_total_vertices);
  yac_remote_point_infos_free(grid.edge_owners, grid.num_total_edges);
  yac_field_data_set_free(grid.field_data);
}
 
void yac_dist_grid_pair_delete(struct dist_grid_pair * grid_pair) {
 
  if (grid_pair == NULL) return;
  yac_mpi_call(MPI_Comm_free(&(grid_pair->comm)), MPI_COMM_WORLD);
  yac_proc_sphere_part_node_delete(grid_pair->proc_sphere_part);
  for (int i = 0; i < 2; ++i) {
    yac_dist_grid_free(grid_pair->dist_grid[i]);
    yac_delete_point_sphere_part_search(grid_pair->vertex_sphere_part[i]);
    yac_bnd_sphere_part_search_delete(grid_pair->cell_sphere_part[i]);
  }
  free(grid_pair);
}
 
static int coord_in_cell(double coord[3], struct dist_grid * grid,
                         size_t cell_idx, struct grid_cell * buffer_cell) {
 
  yac_const_basic_grid_data_get_grid_cell(
    (struct const_basic_grid_data *)grid, cell_idx, buffer_cell);
 
  return
    yac_point_in_cell2(coord, *buffer_cell, grid->cell_bnd_circles[cell_idx]);
}
 
static int coord_in_cell_gc(double coord[3], struct dist_grid * grid,
                         size_t cell_idx, struct grid_cell * buffer_cell) {
 
  yac_const_basic_grid_data_get_grid_cell(
    (struct const_basic_grid_data *)grid, cell_idx, buffer_cell);
  for (size_t i = 0; i < buffer_cell->num_corners; ++i)
    buffer_cell->edge_type[i] = GREAT_CIRCLE_EDGE;
 
  return
    yac_point_in_cell2(coord, *buffer_cell, grid->cell_bnd_circles[cell_idx]);
}
 
void yac_const_basic_grid_data_get_grid_cell(
  struct const_basic_grid_data * grid_data, size_t cell_idx,
  struct grid_cell * buffer_cell) {
 
  size_t num_vertices = (size_t)(grid_data->num_vertices_per_cell[cell_idx]);
 
  struct grid_cell cell = *buffer_cell;
 
  if (cell.array_size < num_vertices) {
    cell.coordinates_xyz =
      xrealloc(cell.coordinates_xyz, num_vertices *
               sizeof(*(cell.coordinates_xyz)));
    cell.edge_type = xrealloc(cell.edge_type, num_vertices *
                                      sizeof(*(cell.edge_type)));
    cell.array_size = num_vertices;
    *buffer_cell = cell;
  }
 
  for (size_t i = 0; i < num_vertices; ++i) {
    size_t vertex_idx =
      grid_data->cell_to_vertex[
        grid_data->cell_to_vertex_offsets[cell_idx] + i];
    cell.coordinates_xyz[i][0] = grid_data->vertex_coordinates[vertex_idx][0];
    cell.coordinates_xyz[i][1] = grid_data->vertex_coordinates[vertex_idx][1];
    cell.coordinates_xyz[i][2] = grid_data->vertex_coordinates[vertex_idx][2];
    size_t edge_idx =
      grid_data->cell_to_edge[grid_data->cell_to_edge_offsets[cell_idx] + i];
    cell.edge_type[i] = grid_data->edge_type[edge_idx];
  }
  buffer_cell->num_corners = num_vertices;
}
 
static struct bnd_sphere_part_search * dist_grid_pair_get_cell_sphere_part(
  struct dist_grid_pair * grid_pair, int grid_id) {
 
  struct bnd_sphere_part_search * search = NULL;
 
  for (int i = 0; (i < 2) && (search == NULL); ++i)
    if (grid_id == grid_pair->grid_ids[i])
      search = grid_pair->cell_sphere_part[i];
  YAC_ASSERT(
    search != NULL,
    "ERROR(yac_dist_grid_pair_get_cell_sphere_part): invalid grid_id")
  return search;
}
 
static void dist_grid_pair_do_point_search_local(
  struct dist_grid_pair * grid_pair, int grid_id,
  coordinate_pointer search_coords, size_t count, size_t * cells,
  int (*coord_in_cell)(
    double coord[3], struct dist_grid * grid, size_t cell_idx,
    struct grid_cell * buffer_cell)) {
 
  struct bnd_sphere_part_search * cell_sphere_part =
    dist_grid_pair_get_cell_sphere_part(grid_pair, grid_id);
  struct dist_grid * grid =
    yac_dist_grid_pair_get_dist_grid(grid_pair, grid_id);
 
  size_t * temp_cells;
  size_t * num_cells_per_coord =
    xmalloc(count * sizeof(*num_cells_per_coord));
 
  // search for all matching source cells
  yac_bnd_sphere_part_search_do_point_search(
    cell_sphere_part, search_coords, count, &temp_cells,
    num_cells_per_coord);
 
  struct grid_cell buffer_cell;
  yac_init_grid_cell(&buffer_cell);
 
  // if we have multiple source cells for a single search coordinate, get the
  // source cell with the lowest global id
  for (size_t i = 0, k = 0; i < count; ++i) {
    size_t curr_num_cells = num_cells_per_coord[i];
    if (curr_num_cells == 0) {
      cells[i] = SIZE_MAX;
    } else if (curr_num_cells == 1) {
      if (coord_in_cell(
            search_coords[i], grid, temp_cells[k], &buffer_cell))
        cells[i] = temp_cells[k];
      else
        cells[i] = SIZE_MAX;
      ++k;
    } else {
      size_t cell_idx = SIZE_MAX;
      yac_int cell_id = XT_INT_MAX;
      for (size_t j = 0; j < curr_num_cells; ++j, ++k) {
        size_t curr_cell_idx = temp_cells[k];
        yac_int curr_cell_id = grid->cell_ids[curr_cell_idx];
        if (!coord_in_cell(
               search_coords[i], grid, curr_cell_idx, &buffer_cell))
          continue;
        if (curr_cell_id < cell_id) {
          cell_idx = curr_cell_idx;
          cell_id = curr_cell_id;
        }
      }
      cells[i] = cell_idx;
    }
  }
 
  yac_free_grid_cell(&buffer_cell);
  free(num_cells_per_coord);
  free(temp_cells);
}
 
static void lookup_single_remote_point_reorder_locally(
  struct dist_grid * grid, enum yac_location location,
  struct single_remote_point_reorder * ids, size_t * count, size_t * idx) {
 
  yac_int * sorted_ids;
  size_t * reorder_idx;
  size_t num_ids;
 
  YAC_ASSERT(
    (location == CELL) || (location == CORNER) || (location == EDGE),
    "ERROR(lookup_single_remote_point_reorder_locally): invalid location")
 
  // get required grid data
  switch (location) {
    default:
    case (CELL):
      sorted_ids = grid->sorted_cell_ids;
      reorder_idx = grid->sorted_cell_reorder_idx;
      num_ids = grid->num_total_cells;
      break;
    case (CORNER):
      sorted_ids = grid->sorted_vertex_ids;
      reorder_idx = grid->sorted_vertex_reorder_idx;
      num_ids = grid->num_total_vertices;
      break;
    case (EDGE):
      sorted_ids = grid->sorted_edge_ids;
      reorder_idx = grid->sorted_edge_reorder_idx;
      num_ids = grid->num_total_edges;
      break;
  };
 
  size_t count_ = *count;
  size_t new_count = 0;
 
  // sort ids by global ids
  qsort(ids, count_, sizeof(*ids),
    compare_single_remote_point_reorder_global_id);
 
  for (size_t i = 0, j = 0; i < count_; ++i) {
    yac_int curr_id = ids[i].data.global_id;
    while ((j < num_ids) && (sorted_ids[j] < curr_id)) ++j;
    if ((j < num_ids) && (sorted_ids[j] == curr_id)) {
      idx[ids[i].reorder_idx] = reorder_idx[j];
    } else {
      if (i != new_count) ids[new_count] = ids[i];
      ++new_count;
    }
  }
 
  *count = new_count;
}
 
static int get_pack_size_field_data(
  struct yac_field_data field_data, MPI_Comm comm) {
 
  int pack_size_field_coord, pack_size_field_mask;
 
  yac_mpi_call(
    MPI_Pack_size(3, MPI_DOUBLE, comm, &pack_size_field_coord), comm);
  pack_size_field_coord *= (int)field_data.coordinates_count;
 
  yac_mpi_call(
    MPI_Pack_size(1, MPI_INT, comm, &pack_size_field_mask), comm);
  pack_size_field_mask *= (int)field_data.masks_count;
 
  return pack_size_field_coord + pack_size_field_mask;
}
 
static int get_pack_size_base_cell(
  struct yac_field_data cell_field_data,
  MPI_Datatype bnd_circle_dt, MPI_Comm comm) {
 
  int pack_size_id,
      pack_size_num_vertices,
      pack_size_bnd_circle;
 
  // id
  yac_mpi_call(MPI_Pack_size(1, Xt_int_dt, comm, &pack_size_id), comm);
  // num_vertices
  yac_mpi_call(MPI_Pack_size(1, MPI_INT, comm, &pack_size_num_vertices), comm);
  // bounding circle
  yac_mpi_call(
    MPI_Pack_size(1, bnd_circle_dt, comm, &pack_size_bnd_circle), comm);
 
  return pack_size_id + pack_size_num_vertices + pack_size_bnd_circle +
         get_pack_size_field_data(cell_field_data, comm);
}
 
static int get_pack_size_base_vertex(
  struct yac_field_data vertex_field_data, MPI_Comm comm) {
 
  int pack_size_id,
      pack_size_vertex_coords;
  // id
  yac_mpi_call(MPI_Pack_size(1, Xt_int_dt, comm, &pack_size_id), comm);
  // vertex coordinates
  yac_mpi_call(
    MPI_Pack_size(3, MPI_DOUBLE, comm, &pack_size_vertex_coords), comm);
 
  return pack_size_id + pack_size_vertex_coords +
         get_pack_size_field_data(vertex_field_data, comm);
}
 
static int get_pack_size_base_edge(
  struct yac_field_data edge_field_data, MPI_Comm comm) {
 
  int pack_size_id,
      pack_size_edge_to_vertex,
      pack_size_edge_type;
  // id
  yac_mpi_call(MPI_Pack_size(1, Xt_int_dt, comm, &pack_size_id), comm);
  // edge type
  yac_mpi_call(MPI_Pack_size(1, MPI_INT, comm, &pack_size_edge_type), comm);
  // edge vertex ids
  yac_mpi_call(
    MPI_Pack_size(2, Xt_int_dt, comm, &pack_size_edge_to_vertex), comm);
 
  return pack_size_id + pack_size_edge_type + pack_size_edge_to_vertex +
         get_pack_size_field_data(edge_field_data, comm);
}
 
static void get_pack_sizes_cell(
  struct dist_grid * grid, uint64_t * pos, size_t count, int * pack_sizes,
  MPI_Datatype bnd_circle_dt, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  int pack_size_base_cell =
    get_pack_size_base_cell(grid->field_data.cell, bnd_circle_dt, comm);
  int pack_size_base_vertex =
    get_pack_size_base_vertex(grid->field_data.vertex, comm);
  int pack_size_base_edge =
    get_pack_size_base_edge(grid->field_data.edge, comm);
 
  for (size_t i = 0; i < count; ++i) {
    size_t idx = (size_t)(pos[i]);
    int num_vertices = grid->num_vertices_per_cell[idx];
    size_t * curr_vertices =
      grid->cell_to_vertex + grid->cell_to_vertex_offsets[idx];
    size_t * curr_edges =
      grid->cell_to_edge + grid->cell_to_edge_offsets[idx];
    int pack_size =
      pack_size_base_cell +
      num_vertices * (pack_size_base_vertex + pack_size_base_edge) +
      yac_remote_point_infos_get_pack_size(
        grid->cell_owners + idx, point_info_dt, comm);
    for (int j = 0; j < num_vertices; ++j) {
      pack_size +=
        yac_remote_point_infos_get_pack_size(
          grid->vertex_owners + curr_vertices[j], point_info_dt, comm) +
        yac_remote_point_infos_get_pack_size(
          grid->edge_owners + curr_edges[j], point_info_dt, comm);
    }
    pack_sizes[i] = pack_size;
  }
}
 
static void get_pack_sizes_vertex(
  struct dist_grid * grid, uint64_t * pos, size_t count, int * pack_sizes,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  int pack_size_base_vertex =
    get_pack_size_base_vertex(grid->field_data.vertex, comm);
  for (size_t i = 0; i < count; ++i)
    pack_sizes[i] =
      pack_size_base_vertex +
      yac_remote_point_infos_get_pack_size(
        grid->vertex_owners + pos[i], point_info_dt, comm);
}
 
static void get_pack_sizes_edge(
  struct dist_grid * grid, uint64_t * pos, size_t count, int * pack_sizes,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  int pack_size_base_vertex =
    get_pack_size_base_vertex(grid->field_data.vertex, comm);
  int pack_size_base_edge =
    get_pack_size_base_edge(grid->field_data.edge, comm);
  for (size_t i = 0; i < count; ++i) {
    size_t * curr_vertices = grid->edge_to_vertex[pos[i]];
    pack_sizes[i] =
      pack_size_base_edge +
      2 * pack_size_base_vertex +
      yac_remote_point_infos_get_pack_size(
        grid->edge_owners + pos[i], point_info_dt, comm) +
      yac_remote_point_infos_get_pack_size(
        grid->vertex_owners + curr_vertices[0], point_info_dt, comm) +
      yac_remote_point_infos_get_pack_size(
        grid->vertex_owners + curr_vertices[1], point_info_dt, comm);
  }
}
 
static void get_pack_sizes(
  struct dist_grid * grid, enum yac_location location, uint64_t * pos,
  size_t count, int * pack_sizes, MPI_Datatype bnd_circle_dt,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  YAC_ASSERT(
    (location == CELL) || (location == CORNER) || (location == EDGE),
    "ERROR(get_pack_sizes): invalid location")
 
  switch(location) {
    default:
    case(CELL):
      get_pack_sizes_cell(
        grid, pos, count, pack_sizes, bnd_circle_dt, point_info_dt, comm);
      break;
    case(CORNER):
      get_pack_sizes_vertex(
        grid, pos, count, pack_sizes, point_info_dt, comm);
      break;
    case(EDGE):
      get_pack_sizes_edge(
        grid, pos, count, pack_sizes, point_info_dt, comm);
      break;
  };
}
 
static void pack_field_data(
  size_t idx, void * buffer, int buffer_size, int * position,
  struct yac_field_data field_data, MPI_Comm comm) {
 
  // coordinates
  for (size_t i = 0; i < field_data.coordinates_count; ++i)
    yac_mpi_call(
      MPI_Pack(field_data.coordinates[i][idx], 3, MPI_DOUBLE, buffer,
               buffer_size, position, comm), comm);
 
  // masks
  for (size_t i = 0; i < field_data.masks_count; ++i)
    yac_mpi_call(
      MPI_Pack(field_data.masks[i] + idx, 1, MPI_INT, buffer,
               buffer_size, position, comm), comm);
}
 
static void pack_grid_data_vertex(
  struct dist_grid * grid, size_t idx, void * buffer, int buffer_size,
  int * position, MPI_Datatype bnd_circle_dt, MPI_Datatype point_info_dt,
  MPI_Comm comm) {
 
  UNUSED(bnd_circle_dt);
 
  // id
  yac_mpi_call(
    MPI_Pack(grid->vertex_ids + idx, 1, Xt_int_dt, buffer, buffer_size,
             position, comm), comm);
  // vertex coordinates
  yac_mpi_call(
    MPI_Pack(&(grid->vertex_coordinates[idx][0]), 3, MPI_DOUBLE, buffer,
             buffer_size, position, comm), comm);
  // vertex owner
  yac_remote_point_infos_pack_single(
    grid->vertex_owners + idx, buffer, buffer_size, position,
    point_info_dt, comm);
  // pack field data
  pack_field_data(
    idx, buffer, buffer_size, position, grid->field_data.vertex, comm);
}
 
static void pack_grid_data_edge_(
  struct dist_grid * grid, size_t idx, void * buffer, int buffer_size,
  int * position, MPI_Datatype bnd_circle_dt, MPI_Datatype point_info_dt,
  MPI_Comm comm) {
 
  UNUSED(bnd_circle_dt);
 
  int edge_type = (int)(grid->edge_type[idx]);
 
  // id
  yac_mpi_call(
    MPI_Pack(grid->edge_ids + idx, 1, Xt_int_dt, buffer, buffer_size, position,
             comm), comm);
  // edge type
  yac_mpi_call(
    MPI_Pack(
      &edge_type, 1, MPI_INT, buffer, buffer_size, position, comm), comm);
  // edge to vertex
  yac_int edge_to_vertex[2] = {
    grid->vertex_ids[grid->edge_to_vertex[idx][0]],
    grid->vertex_ids[grid->edge_to_vertex[idx][1]]};
  yac_mpi_call(
    MPI_Pack(
      edge_to_vertex, 2, Xt_int_dt, buffer, buffer_size, position, comm), comm);
  // edge owner
  yac_remote_point_infos_pack_single(
    grid->edge_owners + idx, buffer, buffer_size, position,
    point_info_dt, comm);
  // pack field data
  pack_field_data(
    idx, buffer, buffer_size, position, grid->field_data.edge, comm);
}
 
static void pack_grid_data_edge(
  struct dist_grid * grid, size_t idx, void * buffer, int buffer_size,
  int * position, MPI_Datatype bnd_circle_dt, MPI_Datatype point_info_dt,
  MPI_Comm comm) {
 
  pack_grid_data_edge_(
    grid, idx, buffer, buffer_size, position,
    bnd_circle_dt, point_info_dt, comm);
 
  // pack edge vertices
  for (int i = 0; i < 2; ++i)
    pack_grid_data_vertex(
      grid, grid->edge_to_vertex[idx][i],
      buffer, buffer_size, position, bnd_circle_dt, point_info_dt, comm);
}
 
static void pack_grid_data_cell(
  struct dist_grid * grid, size_t idx, void * buffer, int buffer_size,
  int * position, MPI_Datatype bnd_circle_dt, MPI_Datatype point_info_dt,
  MPI_Comm comm) {
 
  int num_vertices = grid->num_vertices_per_cell[idx];
 
  // id
  yac_mpi_call(
    MPI_Pack(grid->cell_ids + idx, 1, Xt_int_dt, buffer, buffer_size, position,
             comm), comm);
  // pack field data
  pack_field_data(
    idx, buffer, buffer_size, position, grid->field_data.cell, comm);
  // num_vertices
  yac_mpi_call(
    MPI_Pack(&num_vertices, 1, MPI_INT, buffer,
             buffer_size, position, comm), comm);
  // bounding_circle
  yac_mpi_call(
    MPI_Pack(grid->cell_bnd_circles + idx, 1, bnd_circle_dt, buffer,
             buffer_size, position, comm), comm);
  // cell owner
  yac_remote_point_infos_pack_single(
    grid->cell_owners + idx, buffer, buffer_size, position,
    point_info_dt, comm);
 
  for (int i = 0; i < num_vertices; ++i) {
    pack_grid_data_vertex(
      grid, grid->cell_to_vertex[grid->cell_to_vertex_offsets[idx] + i],
      buffer, buffer_size, position, bnd_circle_dt, point_info_dt, comm);
    pack_grid_data_edge_(
      grid, grid->cell_to_edge[grid->cell_to_edge_offsets[idx] + i],
      buffer, buffer_size, position, bnd_circle_dt, point_info_dt, comm);
  }
}
 
static void pack_grid_data(
  struct dist_grid * grid, enum yac_location location, uint64_t * pos,
  size_t count, void ** pack_data, int * pack_sizes,
  MPI_Datatype bnd_circle_dt, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  get_pack_sizes(grid, location, pos, count, pack_sizes,
                 bnd_circle_dt, point_info_dt, comm);
 
  size_t pack_size = 0;
  for (size_t i = 0; i < count; ++i) pack_size += (size_t)(pack_sizes[i]);
 
  void * pack_data_ = xmalloc(pack_size);
 
  void (*func_pack)(
    struct dist_grid * grid, size_t idx, void * buffer, int buffer_size,
    int * position, MPI_Datatype bnd_circle_dt, MPI_Datatype point_info_dt,
    MPI_Comm comm);
 
  YAC_ASSERT(
    (location == CELL) || (location == CORNER) || (location == EDGE),
    "ERROR(pack_grid_data): invalid location")
 
  switch(location) {
    default:
    case(CELL):
      func_pack = pack_grid_data_cell;
      break;
    case(CORNER):
      func_pack = pack_grid_data_vertex;
      break;
    case(EDGE):
      func_pack = pack_grid_data_edge;
      break;
  };
 
  for (size_t i = 0, offset = 0; i < count; ++i) {
    int position = 0;
    func_pack(grid, pos[i], (char*)pack_data_ + offset, pack_sizes[i],
              &position, bnd_circle_dt, point_info_dt, comm);
    pack_sizes[i] = position;
    offset += (size_t)position;
  }
 
  *pack_data = pack_data_;
}
 
static void unpack_field_data(
  void * buffer, int buffer_size, int * position, size_t idx,
  struct temp_field_data temp_field_data, MPI_Comm comm) {
 
  for (size_t i = 0; i < temp_field_data.coordinates_count; ++i)
    yac_mpi_call(
      MPI_Unpack(
        buffer, buffer_size, position, temp_field_data.coordinates[i][idx],
        3, MPI_DOUBLE, comm), comm);
 
  for (size_t i = 0; i < temp_field_data.masks_count; ++i)
    yac_mpi_call(
      MPI_Unpack(buffer, buffer_size, position,
                 temp_field_data.masks[i] + idx, 1, MPI_INT, comm), comm);
}
 
static void unpack_grid_data_vertex(
  struct global_vertex_reorder * vertex, size_t idx, void * buffer,
  int buffer_size, int * position,
  struct temp_field_data temp_vertex_field_data,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  // id
  yac_mpi_call(
    MPI_Unpack(buffer, buffer_size, position, &(vertex[idx].global_id), 1,
               Xt_int_dt, comm), comm);
  // vertex coordinates
  yac_mpi_call(
    MPI_Unpack(buffer, buffer_size, position, &(vertex[idx].coord[0]), 3,
               MPI_DOUBLE, comm), comm);
  // vertex owners
  yac_remote_point_infos_unpack(
    buffer, buffer_size, position, &(vertex[idx].owners), point_info_dt, comm);
  // unpack field data
  unpack_field_data(
    buffer, buffer_size, position, idx, temp_vertex_field_data, comm);
}
 
static void unpack_grid_data_edge(
  struct global_edge_reorder * edge, size_t idx, void * buffer,
  int buffer_size, int * position,
  struct temp_field_data temp_edge_field_data,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  int edge_type;
 
  // id
  yac_mpi_call(
    MPI_Unpack(buffer, buffer_size, position, &(edge[idx].global_id), 1,
               Xt_int_dt, comm), comm);
  // edge type
  yac_mpi_call(
    MPI_Unpack(buffer, buffer_size, position, &edge_type, 1,
               MPI_INT, comm), comm);
  edge[idx].edge_type = (enum yac_edge_type)edge_type;
  // edge to vertex
  yac_mpi_call(
    MPI_Unpack(buffer, buffer_size, position, edge[idx].edge_to_vertex, 2,
               Xt_int_dt, comm), comm);
  // edge owners
  yac_remote_point_infos_unpack(buffer, buffer_size, position,
                                &(edge[idx].owners), point_info_dt, comm);
  // unpack field data
  unpack_field_data(
    buffer, buffer_size, position, idx, temp_edge_field_data, comm);
}
 
static int compare_global_vertex_reorder_global_id(
  const void * a, const void * b) {
 
  return (((const struct global_vertex_reorder *)a)->global_id >
          ((const struct global_vertex_reorder *)b)->global_id) -
         (((const struct global_vertex_reorder *)a)->global_id <
          ((const struct global_vertex_reorder *)b)->global_id);
}
 
static void add_field_data(
  struct yac_field_data * field_data, struct temp_field_data temp_field_data,
  void * reorder_idx, size_t reorder_idx_size,
  size_t old_count, size_t new_count) {
 
  size_t add_count = new_count - old_count;
 
  for (size_t i = 0; i < temp_field_data.masks_count; ++i) {
    int * temp_mask = temp_field_data.masks[i];
    int * mask =
      ((field_data->masks[i] =
         xrealloc(field_data->masks[i], new_count * sizeof(*mask))));
    for (size_t i = 0, j = old_count; i < add_count; ++i, ++j) {
      size_t idx =
        *(size_t*)((unsigned char*)reorder_idx + i * reorder_idx_size);
      mask[j] = temp_mask[idx];
    }
  }
 
  for (size_t i = 0; i < temp_field_data.coordinates_count; ++i) {
    coordinate_pointer temp_coordinates = temp_field_data.coordinates[i];
    coordinate_pointer coordinates =
      ((field_data->coordinates[i] =
         xrealloc(
          field_data->coordinates[i], new_count * sizeof(*coordinates))));
    for (size_t i = 0, j = old_count; i < add_count; ++i, ++j) {
      size_t idx =
        *(size_t*)((unsigned char*)reorder_idx + i * reorder_idx_size);
      coordinates[j][0] = temp_coordinates[idx][0];
      coordinates[j][1] = temp_coordinates[idx][1];
      coordinates[j][2] = temp_coordinates[idx][2];
    }
  }
}
 
static void yac_remote_point_infos_single_free(
  struct remote_point_infos * point_infos) {
 
  if (point_infos->count > 1) free(point_infos->data.multi);
}
 
static void yac_dist_grid_add_vertices(
  struct dist_grid * grid, struct global_vertex_reorder * vertices,
  size_t count, size_t * idx,
  struct temp_field_data temp_vertex_field_data) {
 
  if (count == 0) return;
 
  // sort vertices global ids
  qsort(vertices, count, sizeof(*vertices),
        compare_global_vertex_reorder_global_id);
 
  yac_int * sorted_vertex_ids = grid->sorted_vertex_ids;
  size_t * sorted_vertex_reorder_idx = grid->sorted_vertex_reorder_idx;
 
  yac_int prev_global_id = vertices[0].global_id - 1;
  size_t prev_idx = 0;
  size_t add_count = 0;
  size_t num_total_vertices = grid->num_total_vertices;
 
  // determine which vertices need to be added to local data
  for (size_t i = 0, j = 0; i < count; ++i) {
 
    yac_int curr_global_id = vertices[i].global_id;
    size_t curr_reorder_idx = vertices[i].reorder_idx;
 
    // if the current global id is a duplicate
    if (prev_global_id == curr_global_id) {
      if (idx != NULL) idx[curr_reorder_idx] = prev_idx;
      yac_remote_point_infos_single_free(&(vertices[i].owners));
      continue;
    }
    prev_global_id = curr_global_id;
 
    // check whether the current global id is already part of the local
    // grid data
    while ((j < num_total_vertices) && (sorted_vertex_ids[j] < curr_global_id))
      ++j;
 
    // if we found a match in the local data
    if ((j < num_total_vertices) && (sorted_vertex_ids[j] == curr_global_id)) {
 
      if (idx != NULL) idx[curr_reorder_idx] = sorted_vertex_reorder_idx[j];
      prev_idx = sorted_vertex_reorder_idx[j];
      yac_remote_point_infos_single_free(&(vertices[i].owners));
 
    // if we need to add the current vertex to the local data
    } else {
 
      if (idx != NULL) idx[curr_reorder_idx] = num_total_vertices + add_count;
      prev_idx = num_total_vertices + add_count;
      if (add_count != i) vertices[add_count] = vertices[i];
      ++add_count;
    }
  }
 
  size_t new_num_total_vertices = num_total_vertices + add_count;
  coordinate_pointer vertex_coordinates =
    xrealloc(grid->vertex_coordinates, new_num_total_vertices *
             sizeof(*vertex_coordinates));
  yac_int * vertex_ids =
    xrealloc(grid->vertex_ids, new_num_total_vertices * sizeof(*vertex_ids));
  int * core_vertex_mask =
    xrealloc(grid->core_vertex_mask, new_num_total_vertices *
             sizeof(*core_vertex_mask));
  struct remote_point_infos * vertex_owners =
    xrealloc(grid->vertex_owners, new_num_total_vertices *
             sizeof(*vertex_owners));
  sorted_vertex_ids =
    xrealloc(
      sorted_vertex_ids, new_num_total_vertices * sizeof(*sorted_vertex_ids));
  sorted_vertex_reorder_idx =
    xrealloc(
      sorted_vertex_reorder_idx, new_num_total_vertices *
      sizeof(*sorted_vertex_reorder_idx));
 
  // add the selected vertices to the local grid data
  for (size_t i = 0, j = num_total_vertices; i < add_count; ++i, ++j) {
 
    vertex_coordinates[j][0] = vertices[i].coord[0];
    vertex_coordinates[j][1] = vertices[i].coord[1];
    vertex_coordinates[j][2] = vertices[i].coord[2];
    vertex_ids[j] = vertices[i].global_id;
    core_vertex_mask[j] = 0;
    vertex_owners[j] = vertices[i].owners;
    sorted_vertex_ids[j] = vertices[i].global_id;
    sorted_vertex_reorder_idx[j] = j;
  }
  // add field data
  add_field_data(
    &(grid->field_data.vertex), temp_vertex_field_data,
    vertices, sizeof(*vertices),
    num_total_vertices, new_num_total_vertices);
  yac_quicksort_index_yac_int_size_t(
    sorted_vertex_ids, new_num_total_vertices, sorted_vertex_reorder_idx);
 
  grid->vertex_coordinates = vertex_coordinates;
  grid->vertex_ids = vertex_ids;
  grid->core_vertex_mask = core_vertex_mask;
  grid->vertex_owners = vertex_owners;
  grid->sorted_vertex_ids = sorted_vertex_ids;
  grid->sorted_vertex_reorder_idx = sorted_vertex_reorder_idx;
  grid->num_total_vertices = new_num_total_vertices;
}
 
static int compare_global_edge_reorder_global_id(
  const void * a, const void * b) {
 
  return (((const struct global_edge_reorder *)a)->global_id >
          ((const struct global_edge_reorder *)b)->global_id) -
         (((const struct global_edge_reorder *)a)->global_id <
          ((const struct global_edge_reorder *)b)->global_id);
}
 
static void yac_dist_grid_add_edges(
  struct dist_grid * grid, struct global_edge_reorder * edges,
  size_t count, size_t * idx, struct temp_field_data temp_edge_field_data) {
 
  if (count == 0) return;
 
  // sort edges global ids
  qsort(edges, count, sizeof(*edges), compare_global_edge_reorder_global_id);
 
  yac_int * sorted_edge_ids = grid->sorted_edge_ids;
  size_t * sorted_edge_reorder_idx = grid->sorted_edge_reorder_idx;
 
  yac_int prev_global_id = edges[0].global_id - 1;
  size_t prev_idx = 0;
  size_t add_count = 0;
  size_t num_total_edges = grid->num_total_edges;
 
  // determine which edges need to be added to local data
  for (size_t i = 0, j = 0; i < count; ++i) {
 
    yac_int curr_global_id = edges[i].global_id;
    size_t curr_reorder_idx = edges[i].reorder_idx;
 
    // if the current global id is a duplicate
    if (prev_global_id == curr_global_id) {
      if (idx != NULL) idx[curr_reorder_idx] = prev_idx;
      yac_remote_point_infos_single_free(&(edges[i].owners));
      continue;
    }
    prev_global_id = curr_global_id;
 
    // check whether the current global id is already part of the local
    // grid data
    while ((j < num_total_edges) && (sorted_edge_ids[j] < curr_global_id)) ++j;
 
    // if we found a match in the local data
    if ((j < num_total_edges) && (sorted_edge_ids[j] == curr_global_id)) {
 
      if (idx != NULL) idx[curr_reorder_idx] = sorted_edge_reorder_idx[j];
      prev_idx = sorted_edge_reorder_idx[j];
      yac_remote_point_infos_single_free(&(edges[i].owners));
 
    // if we need to add the current edge to the local data
    } else {
 
      if (idx != NULL) idx[curr_reorder_idx] = num_total_edges + add_count;
      prev_idx = num_total_edges + add_count;
      if (add_count != i) edges[add_count] = edges[i];
      ++add_count;
    }
  }
 
  size_t new_num_total_edges = num_total_edges + add_count;
  yac_int * edge_ids =
    xrealloc(grid->edge_ids, new_num_total_edges * sizeof(*edge_ids));
  enum yac_edge_type * edge_type =
    xrealloc(grid->edge_type, new_num_total_edges * sizeof(*edge_type));
  size_t_2_pointer edge_to_vertex =
    xrealloc(grid->edge_to_vertex, new_num_total_edges * sizeof(*edge_to_vertex));
  struct remote_point_infos * edge_owners =
    xrealloc(grid->edge_owners, new_num_total_edges * sizeof(*edge_owners));
  int * core_edge_mask =
    xrealloc(grid->core_edge_mask, new_num_total_edges *
             sizeof(*core_edge_mask));
  sorted_edge_ids =
    xrealloc(
      sorted_edge_ids, new_num_total_edges * sizeof(*sorted_edge_ids));
  sorted_edge_reorder_idx =
    xrealloc(
      sorted_edge_reorder_idx, new_num_total_edges *
      sizeof(*sorted_edge_reorder_idx));
 
  yac_int * vertex_ids = xmalloc(2 * add_count * sizeof(*vertex_ids));
  size_t * reorder = xmalloc(2 * add_count * sizeof(*reorder));
 
  // add the selected edges to the local grid data
  for (size_t i = 0, j = num_total_edges; i < add_count; ++i, ++j) {
 
    edge_ids[j] = edges[i].global_id;
    edge_type[j] = edges[i].edge_type;
    core_edge_mask[j] = 0;
    edge_owners[j] = edges[i].owners;
    sorted_edge_ids[j] = edges[i].global_id;
    sorted_edge_reorder_idx[j] = j;
 
    vertex_ids[2 * i + 0] = edges[i].edge_to_vertex[0];
    vertex_ids[2 * i + 1] = edges[i].edge_to_vertex[1];
    reorder[2 * i + 0] = 2 * num_total_edges + 2 * i + 0;
    reorder[2 * i + 1] = 2 * num_total_edges + 2 * i + 1;
  }
  // add field data
  add_field_data(
    &(grid->field_data.edge), temp_edge_field_data, edges, sizeof(*edges),
    num_total_edges, new_num_total_edges);
  yac_quicksort_index_yac_int_size_t(
    sorted_edge_ids, new_num_total_edges, sorted_edge_reorder_idx);
 
  { // determine vertex indices for edge_to_vertex
    yac_quicksort_index_yac_int_size_t(vertex_ids, 2 * add_count, reorder);
    yac_int * sorted_vertex_ids = grid->sorted_vertex_ids;
    size_t * sorted_vertex_reorder_idx = grid->sorted_vertex_reorder_idx;
    size_t total_num_vertices = grid->num_total_vertices;
    size_t * edge_to_vertex_ = (size_t*)&(edge_to_vertex[0][0]);
    // lookup global ids
    for (size_t i = 0, j = 0; i < 2 * add_count; ++i) {
      yac_int curr_id = vertex_ids[i];
      while ((j < total_num_vertices) && (sorted_vertex_ids[j] < curr_id)) ++j;
      YAC_ASSERT(
        (j < total_num_vertices) && (sorted_vertex_ids[j] == curr_id),
        "ERROR(yac_dist_grid_add_edges): vertex id not found")
      edge_to_vertex_[reorder[i]] = sorted_vertex_reorder_idx[j];
    }
  }
 
  free(vertex_ids);
  free(reorder);
 
  grid->edge_ids = edge_ids;
  grid->edge_type = edge_type;
  grid->edge_to_vertex = edge_to_vertex;
  grid->edge_owners = edge_owners;
  grid->core_edge_mask = core_edge_mask;
  grid->sorted_edge_ids = sorted_edge_ids;
  grid->sorted_edge_reorder_idx = sorted_edge_reorder_idx;
  grid->num_total_edges = new_num_total_edges;
}
 
static void yac_dist_grid_add_cells(
  struct dist_grid * grid, yac_int * cell_ids, int * num_vertices_per_cell,
  struct bounding_circle * cell_bnd_circles, size_t count,
  size_t * cell_to_vertex, size_t * cell_to_edge,
  struct remote_point_infos * cell_owners,
  struct temp_field_data temp_cell_field_data) {
 
  if (count == 0) return;
 
  size_t * reorder_idx = xmalloc(count * sizeof(reorder_idx));
  for (size_t i = 0; i < count; ++i) reorder_idx[i] = i;
 
  size_t * prescan = xmalloc(count * sizeof(*prescan));
  for (size_t i = 0, accu = 0; i < count;
       accu += (size_t)(num_vertices_per_cell[i++])) prescan[i] = accu;
 
  // sort cells global ids
  yac_quicksort_index_yac_int_size_t(cell_ids, count, reorder_idx);
 
  yac_int * sorted_cell_ids = grid->sorted_cell_ids;
  size_t * sorted_cell_reorder_idx = grid->sorted_cell_reorder_idx;
 
  yac_int prev_global_id = cell_ids[0] - 1;
  size_t cell_add_count = 0;
  size_t relations_add_count = 0;
  size_t num_total_cells = grid->num_total_cells;
 
  // determine which cells need to be added to local data
  for (size_t i = 0, j = 0; i < count; ++i) {
 
    yac_int curr_global_id = cell_ids[i];
    size_t curr_reorder_idx = reorder_idx[i];
 
    // if the current global id is a duplicate
    if (prev_global_id == curr_global_id) {
      yac_remote_point_infos_single_free(cell_owners + curr_reorder_idx);
      continue;
    }
    prev_global_id = curr_global_id;
 
    // check whether the current global id is already part of the local
    // grid data
    while ((j < num_total_cells) && (sorted_cell_ids[j] < curr_global_id)) ++j;
 
    // if we did not find a match in the local data
    if ((j >= num_total_cells) || (sorted_cell_ids[j] != curr_global_id)) {
 
      if (cell_add_count != i) {
        cell_ids[cell_add_count] = curr_global_id;
        reorder_idx[cell_add_count] = curr_reorder_idx;
      }
      ++cell_add_count;
      relations_add_count += (size_t)(num_vertices_per_cell[i]);
    }
  }
 
  size_t new_num_total_cells = num_total_cells + cell_add_count;
  size_t num_total_relations =
    (num_total_cells > 0)?
      (grid->cell_to_vertex_offsets[num_total_cells-1] +
       (size_t)(grid->num_vertices_per_cell[num_total_cells-1])):0;
  size_t new_num_total_relations = num_total_relations + relations_add_count;
  yac_int * new_cell_ids =
    xrealloc(grid->cell_ids, new_num_total_cells * sizeof(*new_cell_ids));
  int * new_num_vertices_per_cell =
    xrealloc(grid->num_vertices_per_cell, new_num_total_cells *
             sizeof(*new_num_vertices_per_cell));
  size_t * new_cell_to_vertex =
    xrealloc(grid->cell_to_vertex, new_num_total_relations *
             sizeof(*new_cell_to_vertex));
  size_t * cell_to_vertex_offsets =
    xrealloc(grid->cell_to_vertex_offsets, new_num_total_cells *
             sizeof(*cell_to_vertex_offsets));
  size_t * new_cell_to_edge =
    xrealloc(grid->cell_to_edge, new_num_total_relations *
             sizeof(*new_cell_to_edge));
  struct bounding_circle * new_cell_bnd_circles =
    xrealloc(grid->cell_bnd_circles, new_num_total_cells *
             sizeof(*new_cell_bnd_circles));
  int * core_cell_mask =
    xrealloc(grid->core_cell_mask, new_num_total_cells * sizeof(*core_cell_mask));
  struct remote_point_infos * new_cell_owners =
    xrealloc(grid->cell_owners, new_num_total_cells * sizeof(*cell_owners));
  sorted_cell_ids =
    xrealloc(
      sorted_cell_ids, new_num_total_cells * sizeof(*sorted_cell_ids));
  sorted_cell_reorder_idx =
    xrealloc(
      sorted_cell_reorder_idx, new_num_total_cells *
      sizeof(*sorted_cell_reorder_idx));
 
  // add the selected cells to the local grid data
  for (size_t i = 0, j = num_total_cells; i < cell_add_count;
       ++i, ++j) {
 
    size_t curr_reorder_idx = reorder_idx[i];
    int curr_num_vertices = num_vertices_per_cell[curr_reorder_idx];
    size_t curr_relation_idx = prescan[curr_reorder_idx];
 
    new_cell_ids[j] = cell_ids[i];
    new_num_vertices_per_cell[j] = curr_num_vertices;
    cell_to_vertex_offsets[j] = num_total_relations;
    for (int j = 0; j < curr_num_vertices;
         ++j, ++num_total_relations, ++curr_relation_idx) {
      new_cell_to_vertex[num_total_relations] =
        cell_to_vertex[curr_relation_idx];
      new_cell_to_edge[num_total_relations] = cell_to_edge[curr_relation_idx];
    }
    core_cell_mask[j] = 0;
    sorted_cell_ids[j] = cell_ids[i];
    sorted_cell_reorder_idx[j] = j;
    new_cell_bnd_circles[j] = cell_bnd_circles[curr_reorder_idx];
    new_cell_owners[j] = cell_owners[curr_reorder_idx];
  }
  // add field data
  add_field_data(
    &(grid->field_data.cell), temp_cell_field_data,
    reorder_idx, sizeof(*reorder_idx),
    num_total_cells, new_num_total_cells);
  yac_quicksort_index_yac_int_size_t(
    sorted_cell_ids, new_num_total_cells, sorted_cell_reorder_idx);
 
  grid->cell_ids = new_cell_ids;
  grid->num_vertices_per_cell = new_num_vertices_per_cell;
  grid->cell_to_vertex = new_cell_to_vertex;
  grid->cell_to_vertex_offsets = cell_to_vertex_offsets;
  grid->cell_to_edge = new_cell_to_edge;
  grid->cell_to_edge_offsets = cell_to_vertex_offsets;
  grid->core_cell_mask = core_cell_mask;
  grid->cell_owners = new_cell_owners;
  grid->sorted_cell_ids = sorted_cell_ids;
  grid->sorted_cell_reorder_idx = sorted_cell_reorder_idx;
  grid->cell_bnd_circles = new_cell_bnd_circles;
  grid->num_total_cells = new_num_total_cells;
 
  free(prescan);
  free(reorder_idx);
}
 
static void ensure_temp_field_data_sizes(
  struct temp_field_data * temp_field_data, size_t size) {
 
  for (size_t i = 0; i < temp_field_data->masks_count; ++i)
    ENSURE_ARRAY_SIZE(
      temp_field_data->masks[i], temp_field_data->masks_array_sizes[i], size);
  for (size_t i = 0; i < temp_field_data->coordinates_count; ++i)
    ENSURE_ARRAY_SIZE(
      temp_field_data->coordinates[i],
      temp_field_data->coordinates_array_sizes[i], size);
}
 
static struct temp_field_data temp_field_data_init(
  struct yac_field_data field_data, size_t count) {
 
  struct temp_field_data temp_field_data;
 
  temp_field_data.masks =
    xmalloc(field_data.masks_count * sizeof(*temp_field_data.masks));
  temp_field_data.masks_array_sizes =
    xmalloc(
      field_data.masks_count * sizeof(*temp_field_data.masks_array_sizes));
  temp_field_data.masks_count = field_data.masks_count;
  for (size_t i = 0; i < field_data.masks_count; ++i) {
    temp_field_data.masks[i] =
      xmalloc(count * sizeof(**temp_field_data.masks));
    temp_field_data.masks_array_sizes[i] = count;
  }
 
  temp_field_data.coordinates =
    xmalloc(
      field_data.coordinates_count * sizeof(*temp_field_data.coordinates));
  temp_field_data.coordinates_array_sizes =
    xmalloc(
      field_data.coordinates_count *
      sizeof(*temp_field_data.coordinates_array_sizes));
  temp_field_data.coordinates_count = field_data.coordinates_count;
  for (size_t i = 0; i < field_data.coordinates_count; ++i) {
    temp_field_data.coordinates[i] =
      xmalloc(count * sizeof(**temp_field_data.coordinates));
    temp_field_data.coordinates_array_sizes[i] = count;
  }
 
  return temp_field_data;
}
 
static void temp_field_data_free(struct temp_field_data temp_field_data) {
 
  for (size_t i = 0; i < temp_field_data.masks_count; ++i)
    free(temp_field_data.masks[i]);
  free(temp_field_data.masks);
  free(temp_field_data.masks_array_sizes);
  for (size_t i = 0; i < temp_field_data.coordinates_count; ++i)
    free(temp_field_data.coordinates[i]);
  free(temp_field_data.coordinates);
  free(temp_field_data.coordinates_array_sizes);
}
 
static void unpack_grid_data_cells(
  struct dist_grid * grid, size_t count, void * buffer, int buffer_size,
  MPI_Datatype bnd_circle_dt, MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  yac_int * cell_ids = xmalloc(count * sizeof(*cell_ids));
  int * num_vertices_per_cell = xmalloc(count * sizeof(*num_vertices_per_cell));
  struct bounding_circle * cell_bnd_circles =
    xmalloc(count * sizeof(*cell_bnd_circles));
  struct remote_point_infos * cell_owners =
    xmalloc(count * sizeof(*cell_owners));
 
  struct global_vertex_reorder * vertices = NULL;
  size_t vertices_array_size = 0;
  size_t total_num_vertices = 0;
 
  struct global_edge_reorder * edges = NULL;
  size_t edges_array_size = 0;
 
  struct temp_field_data temp_cell_field_data =
    temp_field_data_init(grid->field_data.cell, count);
  struct temp_field_data temp_vertex_field_data =
    temp_field_data_init(grid->field_data.vertex, 3 * count);
  struct temp_field_data temp_edge_field_data =
    temp_field_data_init(grid->field_data.edge, 3 * count);
 
  for (size_t i = 0, buffer_offset = 0; i < count; ++i) {
 
    int position = 0;
    void * curr_buffer = (char*)buffer + buffer_offset;
    int num_vertices;
 
    // cell id
    yac_mpi_call(
      MPI_Unpack(curr_buffer, buffer_size, &position, cell_ids + i, 1,
                 Xt_int_dt, comm), comm);
    // unpack field data
    unpack_field_data(
      curr_buffer, buffer_size, &position, i, temp_cell_field_data, comm);
    // num vertices
    yac_mpi_call(
      MPI_Unpack(curr_buffer, buffer_size, &position, &num_vertices, 1,
                 MPI_INT, comm), comm);
    // bounding circle
    yac_mpi_call(
      MPI_Unpack(curr_buffer, buffer_size, &position, cell_bnd_circles + i, 1,
                 bnd_circle_dt, comm), comm);
    // cell owners
    yac_remote_point_infos_unpack(
      curr_buffer, buffer_size, &position, cell_owners + i,
      point_info_dt, comm);
 
    num_vertices_per_cell[i] = num_vertices;
 
    ENSURE_ARRAY_SIZE(
      vertices, vertices_array_size, total_num_vertices + (size_t)num_vertices);
    ENSURE_ARRAY_SIZE(
      edges, edges_array_size, total_num_vertices + (size_t)num_vertices);
    ensure_temp_field_data_sizes(
      &temp_vertex_field_data, total_num_vertices + (size_t)num_vertices);
    ensure_temp_field_data_sizes(
      &temp_edge_field_data, total_num_vertices + (size_t)num_vertices);
 
    for (int j = 0; j < num_vertices; ++j, ++total_num_vertices) {
      unpack_grid_data_vertex(
        vertices, total_num_vertices, curr_buffer, buffer_size, &position,
        temp_vertex_field_data, point_info_dt, comm);
      unpack_grid_data_edge(
        edges, total_num_vertices, curr_buffer, buffer_size, &position,
        temp_edge_field_data, point_info_dt, comm);
      vertices[total_num_vertices].reorder_idx = total_num_vertices;
      edges[total_num_vertices].reorder_idx = total_num_vertices;
    }
 
    buffer_offset += (size_t)position;
    buffer_size -= position;
  }
 
  size_t * cell_to_vertex = xmalloc(total_num_vertices * sizeof(*cell_to_vertex));
  size_t * cell_to_edge = xmalloc(total_num_vertices * sizeof(*cell_to_edge));
 
  yac_dist_grid_add_vertices(
    grid, vertices, total_num_vertices, cell_to_vertex,
    temp_vertex_field_data);
  yac_dist_grid_add_edges(
    grid, edges, total_num_vertices, cell_to_edge,
    temp_edge_field_data);
  yac_dist_grid_add_cells(
    grid, cell_ids, num_vertices_per_cell, cell_bnd_circles, count,
    cell_to_vertex, cell_to_edge, cell_owners, temp_cell_field_data);
 
  temp_field_data_free(temp_cell_field_data);
  temp_field_data_free(temp_vertex_field_data);
  temp_field_data_free(temp_edge_field_data);
  free(cell_to_edge);
  free(cell_to_vertex);
  free(vertices);
  free(edges);
  free(cell_owners);
  free(cell_bnd_circles);
  free(num_vertices_per_cell);
  free(cell_ids);
}
 
static void unpack_grid_data_vertices(
  struct dist_grid * grid, size_t count, void * buffer, int buffer_size,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  struct global_vertex_reorder * vertices = xmalloc(count * sizeof(*vertices));
 
  struct temp_field_data temp_vertex_field_data =
    temp_field_data_init(grid->field_data.vertex, count);
 
  for (size_t i = 0, buffer_offset = 0; i < count; ++i) {
 
    int position = 0;
    void * curr_buffer = (char*)buffer + buffer_offset;
 
    unpack_grid_data_vertex(
      vertices, i, curr_buffer, buffer_size, &position,
      temp_vertex_field_data, point_info_dt, comm);
    vertices[i].reorder_idx = i;
 
    buffer_offset += (size_t)position;
    buffer_size -= position;
  }
 
  yac_dist_grid_add_vertices(
    grid, vertices, count, NULL, temp_vertex_field_data);
 
  temp_field_data_free(temp_vertex_field_data);
 
  free(vertices);
}
 
static void unpack_grid_data_edges(
  struct dist_grid * grid, size_t count, void * buffer, int buffer_size,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  struct global_edge_reorder * edges = xmalloc(count * sizeof(*edges));
  struct global_vertex_reorder * vertices =
    xmalloc(2 * count * sizeof(*vertices));
 
  struct temp_field_data temp_edge_field_data =
    temp_field_data_init(grid->field_data.edge, count);
  struct temp_field_data temp_vertex_field_data =
    temp_field_data_init(grid->field_data.vertex, 2 * count);
 
  for (size_t i = 0, buffer_offset = 0; i < count; ++i) {
 
    int position = 0;
    void * curr_buffer = (char*)buffer + buffer_offset;
 
    unpack_grid_data_edge(
      edges, i, curr_buffer, buffer_size, &position,
      temp_edge_field_data, point_info_dt, comm);
    edges[i].reorder_idx = i;
 
    for (size_t j = 0; j < 2; ++j)
      unpack_grid_data_vertex(
        vertices, 2 * i + j, curr_buffer, buffer_size, &position,
        temp_vertex_field_data, point_info_dt, comm);
 
    buffer_offset += (size_t)position;
    buffer_size -= position;
  }
 
  yac_dist_grid_add_vertices(
    grid, vertices, 2 * count, NULL, temp_vertex_field_data);
  yac_dist_grid_add_edges(
    grid, edges, count, NULL, temp_edge_field_data);
 
  temp_field_data_free(temp_vertex_field_data);
  temp_field_data_free(temp_edge_field_data);
 
  free(vertices);
  free(edges);
}
 
static void unpack_grid_data(
  struct dist_grid * grid, enum yac_location location, size_t count,
  void * buffer, int buffer_size, MPI_Datatype bnd_circle_dt,
  MPI_Datatype point_info_dt, MPI_Comm comm) {
 
  YAC_ASSERT(
    (location == CELL) || (location == CORNER) || (location == EDGE),
    "ERROR(unpack_grid_data): invalid location")
 
  switch(location) {
    default:
    case(CELL):
      unpack_grid_data_cells(
        grid, count, buffer, buffer_size, bnd_circle_dt, point_info_dt, comm);
      break;
    case(CORNER):
      unpack_grid_data_vertices(
        grid, count, buffer, buffer_size, point_info_dt, comm);
      break;
    case(EDGE):
      unpack_grid_data_edges(
        grid, count, buffer, buffer_size, point_info_dt, comm);
      break;
  };
}
 
static int compare_single_remote_point_reorder_owner(
  const void * a, const void * b) {
 
  return ((const struct single_remote_point_reorder *)a)->data.data.rank -
         ((const struct single_remote_point_reorder *)b)->data.data.rank;
}
 
static void yac_dist_grid_single_remote_point_to_local(
  struct dist_grid * grid, struct single_remote_point * ids, size_t count,
  enum yac_location location, size_t * idx) {
 
  MPI_Comm comm = grid->comm;
  int comm_rank, comm_size;
  yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  size_t remote_count = 0;
 
  for (size_t i = 0; i < count; ++i) {
    if (ids[i].global_id == XT_INT_MAX) idx[i] = SIZE_MAX;
    else if (ids[i].data.rank != comm_rank) ++remote_count;
    else idx[i] = ids[i].data.orig_pos;
  }
 
  struct single_remote_point_reorder * missing_ids =
    xmalloc(remote_count * sizeof(*missing_ids));
 
  for (size_t i = 0, j = 0; i < count; ++i) {
    if ((ids[i].data.rank != comm_rank) &&
        (ids[i].global_id != XT_INT_MAX)) {
      missing_ids[j].data = ids[i];
      missing_ids[j].reorder_idx = i;
      ++j;
    }
  }
 
  // check whether we already have some of the missing ids locally
  lookup_single_remote_point_reorder_locally(
    grid, location, missing_ids, &remote_count, idx);
 
  // sort data by owner
  qsort(missing_ids, remote_count, sizeof(*missing_ids),
        compare_single_remote_point_reorder_owner);
 
  int * sendcounts, * recvcounts, * sdispls, * rdispls;
  yac_get_comm_buffers(
    1, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
 
  for (size_t i = 0; i < remote_count; ++i)
    sendcounts[missing_ids[i].data.data.rank]++;
 
  yac_generate_alltoallv_args(
    1, sendcounts, recvcounts, sdispls, rdispls, comm);
 
  size_t recv_count =
    (size_t)(rdispls[comm_size-1]) + (size_t)(recvcounts[comm_size-1]);
 
  uint64_t * uint64_t_buffer =
    xmalloc((remote_count + recv_count) * sizeof(*uint64_t_buffer));
  uint64_t * orig_pos_send_buffer = uint64_t_buffer;
  uint64_t * orig_pos_recv_buffer = uint64_t_buffer + remote_count;
 
  // pack send buffer
  for (size_t i = 0; i < remote_count; ++i) {
    int rank = missing_ids[i].data.data.rank;
    if (rank != comm_rank)
      orig_pos_send_buffer[sdispls[rank+1]++] =
        (uint64_t)(missing_ids[i].data.data.orig_pos);
  }
 
  // redistribute ids
  yac_alltoallv_uint64_p2p(
    orig_pos_send_buffer, sendcounts, sdispls,
    orig_pos_recv_buffer, recvcounts, rdispls, comm);
 
  MPI_Datatype bnd_circle_dt = yac_get_bounding_circle_mpi_datatype(comm);
  yac_mpi_call(MPI_Type_commit(&bnd_circle_dt), comm);
  MPI_Datatype point_info_dt = yac_get_remote_point_info_mpi_datatype(comm);
  yac_mpi_call(MPI_Type_commit(&point_info_dt), comm);
 
  void * packed_send_data = NULL;
  int * pack_sizes = xmalloc(recv_count * sizeof(*pack_sizes));
 
  // pack all requested grid data
  pack_grid_data(
    grid, location, orig_pos_recv_buffer, recv_count,
    &packed_send_data, pack_sizes,
    bnd_circle_dt, point_info_dt, comm);
  free(uint64_t_buffer);
 
  memset(sendcounts, 0, (size_t)comm_size * sizeof(*sendcounts));
  for (int i = 0, k = 0; i < comm_size; ++i)
    for (int j = 0; j < recvcounts[i]; ++j, ++k)
      sendcounts[i] += pack_sizes[k];
 
  free(pack_sizes);
 
  yac_generate_alltoallv_args(
    1, sendcounts, recvcounts, sdispls, rdispls, comm);
 
  recv_count =
    (size_t)(rdispls[comm_size-1]) + (size_t)(recvcounts[comm_size-1]);
 
  void * packed_recv_data = xmalloc(recv_count);
 
  // redistribute packed grid data
  yac_alltoallv_packed_p2p(
    packed_send_data, sendcounts, sdispls+1,
    packed_recv_data, recvcounts, rdispls, comm);
 
  // unpack requested grid data
  unpack_grid_data(
    grid, location, remote_count, packed_recv_data, (int)recv_count,
    bnd_circle_dt, point_info_dt, comm);
 
  yac_mpi_call(MPI_Type_free(&point_info_dt), comm);
  yac_mpi_call(MPI_Type_free(&bnd_circle_dt), comm);
 
  // get the local ids for the remaining missing ids
  lookup_single_remote_point_reorder_locally(
    grid, location, missing_ids, &remote_count, idx);
 
  free(missing_ids);
  free(packed_recv_data);
  free(packed_send_data);
  yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
}
 
void yac_dist_grid_pair_do_point_search_(
  struct dist_grid_pair * grid_pair, int grid_id,
  coordinate_pointer search_coords, size_t count, size_t * cells,
  int (*coord_in_cell)(
    double coord[3], struct dist_grid * grid, size_t cell_idx,
    struct grid_cell * buffer_cell)) {
 
  MPI_Comm comm = grid_pair->comm;
  int comm_rank, comm_size;
  yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  int * ranks = xmalloc(count * sizeof(ranks));
 
  // search for every point in the process tree
  yac_proc_sphere_part_do_point_search(
    grid_pair->proc_sphere_part, search_coords, count, ranks);
 
  int * sendcounts, * recvcounts, * sdispls, * rdispls;
  yac_get_comm_buffers(
    1, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
  for (size_t i = 0; i < count; ++i) sendcounts[ranks[i]]++;
 
  size_t local_count = sendcounts[comm_rank];
  sendcounts[comm_rank] = 0;
 
  yac_generate_alltoallv_args(
    1, sendcounts, recvcounts, sdispls, rdispls, comm);
 
  size_t remote_count = sdispls[comm_size] + sendcounts[comm_size-1];
  size_t request_count = rdispls[comm_size-1] + recvcounts[comm_size-1];
 
  coordinate_pointer coord_buffer =
    xmalloc((remote_count + request_count + local_count) *
            sizeof(*coord_buffer));
  coordinate_pointer coord_send_buffer = coord_buffer + 0;
  coordinate_pointer coord_recv_buffer = coord_buffer + remote_count;
  coordinate_pointer coord_local_buffer =
    coord_buffer + remote_count + request_count;
 
  // pack search coordinates
  for (size_t i = 0, k = 0; i < count; ++i) {
    if (ranks[i] == comm_rank) {
      coord_local_buffer[k][0] = search_coords[i][0];
      coord_local_buffer[k][1] = search_coords[i][1];
      coord_local_buffer[k][2] = search_coords[i][2];
      ++k;
    } else {
      int displ = sdispls[ranks[i]+1]++;
      coord_send_buffer[displ][0] = search_coords[i][0];
      coord_send_buffer[displ][1] = search_coords[i][1];
      coord_send_buffer[displ][2] = search_coords[i][2];
    }
  }
 
  MPI_Datatype dt_coord;
  yac_mpi_call(MPI_Type_contiguous(3, MPI_DOUBLE, &dt_coord), comm);
  yac_mpi_call(MPI_Type_commit(&dt_coord), comm);
 
  // redistribute search coordinates
  yac_alltoallv_p2p(
    coord_send_buffer, sendcounts, sdispls,
    coord_recv_buffer, recvcounts, rdispls,
    sizeof(*coord_send_buffer), dt_coord, comm);
 
  yac_mpi_call(MPI_Type_free(&dt_coord), comm);
 
  size_t * local_cells =
    xmalloc((request_count + local_count) * sizeof(*local_cells));
 
  // do local search
  dist_grid_pair_do_point_search_local(
    grid_pair, grid_id, coord_recv_buffer, request_count + local_count,
    local_cells, coord_in_cell);
 
  struct single_remote_point * single_remote_point_buffer =
    xmalloc((remote_count + request_count) *
            sizeof(*single_remote_point_buffer));
  struct single_remote_point * id_send_buffer = single_remote_point_buffer;
  struct single_remote_point * id_recv_buffer = single_remote_point_buffer +
                                                request_count;
 
  struct dist_grid * grid =
    yac_dist_grid_pair_get_dist_grid(grid_pair, grid_id);
 
  // pack global ids of found source cells
  for (size_t i = 0; i < request_count; ++i) {
    size_t cell_idx = local_cells[i];
    id_send_buffer[i].data.rank = comm_rank;
    if (cell_idx != SIZE_MAX) {
      id_send_buffer[i].global_id = grid->cell_ids[cell_idx];
      id_send_buffer[i].data.orig_pos = (uint64_t)cell_idx;
    } else {
      id_send_buffer[i].global_id = XT_INT_MAX;
      id_send_buffer[i].data.orig_pos = UINT64_MAX;
    }
  }
 
  MPI_Datatype single_remote_point_dt =
    yac_get_single_remote_point_mpi_datatype(comm);
 
  // redistribute results (global ids of found source cells)
  yac_alltoallv_p2p(
    id_send_buffer, recvcounts, rdispls, id_recv_buffer, sendcounts, sdispls,
    sizeof(*id_send_buffer), single_remote_point_dt, comm);
 
  yac_mpi_call(MPI_Type_free(&single_remote_point_dt), comm);
 
  size_t * new_local_cells =
    xmalloc(remote_count * sizeof(*new_local_cells));
 
  // convert all remote ids to local ones, extend local dist_grid data,
  // if necessary
  yac_dist_grid_single_remote_point_to_local(
    grid, id_recv_buffer, remote_count, CELL, new_local_cells);
 
  // extract results from local and remote search
  for (size_t i = 0, k = 0; i < count; ++i) {
    if (ranks[i] == comm_rank) {
      cells[i] = local_cells[request_count + k];
      ++k;
    } else {
      int displ = sdispls[ranks[i]]++;
      cells[i] = new_local_cells[displ];
    }
  }
 
  free(new_local_cells);
  free(single_remote_point_buffer);
  free(local_cells);
  free(coord_buffer);
  yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
  free(ranks);
}
 
void yac_dist_grid_pair_do_point_search(
  struct dist_grid_pair * grid_pair, int grid_id,
  coordinate_pointer search_coords, size_t count, size_t * cells) {
 
  yac_dist_grid_pair_do_point_search_(
    grid_pair, grid_id, search_coords, count, cells, coord_in_cell);
}
 
void yac_dist_grid_pair_do_point_search_gc(
  struct dist_grid_pair * grid_pair, int grid_id,
  coordinate_pointer search_coords, size_t count, size_t * cells) {
 
  yac_dist_grid_pair_do_point_search_(
    grid_pair, grid_id, search_coords, count, cells, coord_in_cell_gc);
}
 
static struct point_sphere_part_search * yac_dist_grid_get_field_sphere_part(
  struct dist_grid * grid, struct interp_field field) {
 
  coordinate_pointer field_coords =
    yac_dist_grid_get_field_coords(grid, field);
  yac_int const * global_ids =
    yac_dist_grid_get_global_ids(grid, field.location);
  int const * mask = yac_dist_grid_get_field_mask(grid, field);
  size_t total_count = yac_dist_grid_get_count(grid, field.location);
 
  if (mask == NULL)
    return
      yac_point_sphere_part_search_new(
        total_count, field_coords, global_ids);
  else
    return
      yac_point_sphere_part_search_mask_new(
        total_count, field_coords, global_ids, mask);
}
 
static void yac_dist_grid_get_n_unmasked_local_points(
  struct dist_grid * grid, struct interp_field field, int comm_rank,
  size_t n, struct single_remote_point * points) {
 
  size_t total_count = yac_dist_grid_get_count(grid, field.location);
  int const * mask = yac_dist_grid_get_field_mask(grid, field);
  int const * core_mask =
    yac_dist_grid_get_core_mask(grid, field.location);
  yac_int const * global_ids =
    yac_dist_grid_get_global_ids(grid, field.location);
 
  if (mask == NULL) {
 
    for (size_t i = 0, j = 0; i < total_count; ++i) {
      if (core_mask[i]) {
        points[j].global_id = global_ids[i];
        points[j].data.rank = comm_rank;
        points[j].data.orig_pos = i;
        if (n == ++j) return;
      }
    }
 
  } else {
 
    for (size_t i = 0, j = 0; i < total_count; ++i) {
      if (core_mask[i] && mask[i]) {
        points[j].global_id = global_ids[i];
        points[j].data.rank = comm_rank;
        points[j].data.orig_pos = i;
        if (n == ++j) return;
      }
    }
  }
}
 
static MPI_Datatype yac_get_nnn_search_result_mpi_datatype(
  size_t n, MPI_Comm comm) {
 
  struct nnn_search_result dummy;
  MPI_Datatype single_remote_point_dt =
                 yac_get_single_remote_point_mpi_datatype(comm),
               single_nnn_search_result_dt_,
               nnn_search_result_dt;
  int array_of_blocklengths[] = {1, 1};
  const MPI_Aint array_of_displacements[] =
    {(MPI_Aint)((void*)&(dummy.point_info) - (void*)&dummy),
     (MPI_Aint)((void*)&(dummy.cos_angle) - (void*)&dummy)};
  const MPI_Datatype array_of_types[] = {single_remote_point_dt, MPI_DOUBLE};
  yac_mpi_call(
    MPI_Type_create_struct(
      2, array_of_blocklengths, array_of_displacements,
      array_of_types, &single_nnn_search_result_dt_), comm);
  MPI_Datatype single_nnn_search_result_dt =
    yac_create_resized(single_nnn_search_result_dt_, sizeof(dummy), comm);
  yac_mpi_call(
    MPI_Type_contiguous(
      (int)n, single_nnn_search_result_dt, &nnn_search_result_dt), comm);
  yac_mpi_call(MPI_Type_commit(&nnn_search_result_dt), comm);
  yac_mpi_call(MPI_Type_free(&single_nnn_search_result_dt), comm);
  yac_mpi_call(MPI_Type_free(&single_remote_point_dt), comm);
  return nnn_search_result_dt;
}
 
static int compare_nnn_search_results_cos_angle(
  void const * a, void const * b) {
 
  struct nnn_search_result const * result_a =
    (struct nnn_search_result const *)a;
  struct nnn_search_result const * result_b =
    (struct nnn_search_result const *)b;
 
  int ret = (result_a->cos_angle < result_b->cos_angle) -
            (result_a->cos_angle > result_b->cos_angle);
  if (ret) return ret;
  return (result_a->point_info.global_id > result_b->point_info.global_id) -
         (result_a->point_info.global_id < result_b->point_info.global_id);
}
 
static void insert_nnn_result_points(
  struct nnn_search_result * array, size_t array_size,
  struct nnn_search_result * insert, size_t insert_size) {
 
  if (insert_size == 0) return;
 
  // sort arrays first by cos_angle and second by global id
  qsort(array, array_size, sizeof(*array),
        compare_nnn_search_results_cos_angle);
  qsort(insert, insert_size, sizeof(*insert),
        compare_nnn_search_results_cos_angle);
 
  struct nnn_search_result temp_results[array_size];
  memcpy(temp_results, array, array_size * sizeof(*array));
 
  for (size_t i = 0, j = 0, k = 0; k < array_size; ++k) {
 
    int cmp =
      (j < insert_size)?
        compare_nnn_search_results_cos_angle(
          &temp_results[i], insert + j):(-1);
 
    if (cmp <= 0) {
      if (k != i) array[k] = temp_results[i];
      ++i;
    } else {
      array[k] = insert[j];
    }
 
    j += cmp >= 0;
  }
}
 
void yac_dist_grid_pair_do_nnn_search(
  struct dist_grid_pair * grid_pair, int grid_id,
  coordinate_pointer search_coords, size_t count, size_t * local_ids,
  size_t n, struct interp_field field) {
 
  MPI_Comm comm = grid_pair->comm;
  int comm_rank, comm_size;
  yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  struct dist_grid * grid =
    yac_dist_grid_pair_get_dist_grid(grid_pair, grid_id);
 
  uint64_t unmasked_local_count =
    yac_dist_grid_get_unmasked_local_count(grid, field);
 
  uint64_t * unmasked_local_counts =
    xmalloc((size_t)comm_size * sizeof(*unmasked_local_counts));
 
  // exchange number of local source points and target points
  yac_mpi_call(
    MPI_Allgather(
      &unmasked_local_count, 1, MPI_UINT64_T,
      unmasked_local_counts, 1, MPI_UINT64_T, comm), comm);
 
  // check whether there is a rank with too few source points
  int flag = 0;
  for (int i = 0; i < comm_size; ++i)
    flag |= unmasked_local_counts[i] < (uint64_t)n;
 
  // if ranks with insufficient number of local source points
  if (flag) {
 
    uint64_t global_num_unmasked_count = 0;
    for (int i = 0; i < comm_size; ++i)
      global_num_unmasked_count += unmasked_local_counts[i];
 
    YAC_ASSERT(
      (size_t)global_num_unmasked_count >= n,
      "ERROR(yac_dist_grid_pair_do_nnn_search): insufficient number of "
      "unmasked points")
 
    int * sendcounts, * recvcounts, * sdispls, * rdispls;
    yac_get_comm_buffers(
      1, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
    memset(recvcounts, 0, (size_t)comm_size * sizeof(*recvcounts));
 
    // get ranks of processes that have additional data or require data in
    // order to have enough source points to do a initial nnn search
    yac_proc_sphere_part_get_neigh_ranks(
      grid_pair->proc_sphere_part, unmasked_local_counts, (uint64_t)n,
      sendcounts, recvcounts, comm_rank, comm_size);
 
    size_t local_send_count = (size_t)(MIN(unmasked_local_count, n));
 
    for (int i = 0, raccu = 0; i < comm_size; ++i) {
      sdispls[i] = 0;
      rdispls[i] = raccu;
      sendcounts[i] *= (int)local_send_count;
      raccu += (recvcounts[i] *= (int)(MIN(unmasked_local_counts[i], n)));
    }
 
    size_t recv_count =
      (size_t)(recvcounts[comm_size-1]) + (size_t)(rdispls[comm_size-1]);
 
    struct single_remote_point * single_remote_point_buffer =
      xmalloc(
        (local_send_count + recv_count) * sizeof(*single_remote_point_buffer));
    struct single_remote_point * local_send_ids = single_remote_point_buffer;
    struct single_remote_point * recv_ids =
      single_remote_point_buffer + local_send_count;
 
    // get local source points that can be sent to other processes
    yac_dist_grid_get_n_unmasked_local_points(
      grid, field, comm_rank, local_send_count, local_send_ids);
 
    MPI_Datatype single_remote_point_dt =
      yac_get_single_remote_point_mpi_datatype(comm);
 
    // exchange source points (integrate points into local data)
    yac_alltoallv_p2p(
      local_send_ids, sendcounts, sdispls, recv_ids, recvcounts, rdispls,
      sizeof(*local_send_ids), single_remote_point_dt, comm);
 
    yac_mpi_call(MPI_Type_free(&single_remote_point_dt), comm);
 
    size_t * dummy = xmalloc(recv_count * sizeof(*dummy));
 
    // convert all remote ids to local ones, extend local dist_grid data,
    // if necessary
    yac_dist_grid_single_remote_point_to_local(
      grid, recv_ids, recv_count, field.location, dummy);
 
    free(dummy);
    free(single_remote_point_buffer);
    yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
  }
  free(unmasked_local_counts);
 
  struct point_sphere_part_search * sphere_part =
    yac_dist_grid_get_field_sphere_part(grid, field);
 
  double * cos_angles = NULL;
  size_t cos_angles_array_size = 0;
  coordinate_pointer result_coords = NULL;
  size_t result_coords_array_size = 0;
  size_t * result_points = NULL;
  size_t result_points_array_size = 0;
  size_t * num_results = xcalloc(count, sizeof(*num_results));
 
  // do local search
  yac_point_sphere_part_search_NNN(
    sphere_part, count, search_coords, n, &cos_angles, &cos_angles_array_size,
    &result_coords, &result_coords_array_size, &result_points,
    &result_points_array_size, num_results);
 
  size_t max_num_results = 0;
  for (size_t i = 0; i < count; ++i)
    if (max_num_results < num_results[i]) max_num_results = num_results[i];
 
  // extract results (only the local core points)
  struct nnn_search_result * results = xmalloc(n * count * sizeof(*results));
  yac_int const * global_ids =
    yac_dist_grid_get_global_ids(grid, field.location);
 
  struct nnn_search_result * temp_results =
    xmalloc(max_num_results * sizeof(*temp_results));
 
  for (size_t i = 0, k = 0; i < count; ++i) {
 
    struct nnn_search_result * curr_results = results + i * n;
    for (size_t j = 0; j < n; ++j, ++k) {
      size_t curr_local_id = result_points[k];
      curr_results[j].point_info.global_id = global_ids[curr_local_id];
      curr_results[j].point_info.data.rank = comm_rank;
      curr_results[j].point_info.data.orig_pos = (uint64_t)curr_local_id;
      curr_results[j].cos_angle = cos_angles[k];
    }
 
    for (size_t j = n; j < num_results[i]; ++j, ++k) {
      size_t curr_local_id = result_points[k];
      temp_results[j-n].point_info.global_id = global_ids[curr_local_id];
      temp_results[j-n].point_info.data.rank = comm_rank;
      temp_results[j-n].point_info.data.orig_pos = (uint64_t)curr_local_id;
      temp_results[j-n].cos_angle = cos_angles[k];
    }
    insert_nnn_result_points(curr_results, n, temp_results, num_results[i] - n);
  }
 
  coordinate_pointer field_coords =
    yac_dist_grid_get_field_coords(grid, field);
  int * request_ranks = NULL;
  size_t request_ranks_array_size = 0;
  size_t num_request_ranks = 0;
  int * num_requests = xmalloc(count * sizeof(*num_requests));
 
  // for all search points
  for (size_t i = 0; i < count; ++i) {
 
    // generate bounding circles for all search points
    struct bounding_circle result_bnd_circle;
    memcpy(result_bnd_circle.base_vector, search_coords[i],
           3 * sizeof(search_coords[0][0]));
    result_bnd_circle.inc_angle =
      get_vector_angle_2(
        search_coords[i],
        field_coords[results[(i + 1) * n - 1].point_info.data.orig_pos]);
 
    ENSURE_ARRAY_SIZE(request_ranks, request_ranks_array_size,
                      num_request_ranks + (size_t)comm_size);
 
    // search for processes that might be able to contribute to the search
    // results
    int * curr_request_ranks = request_ranks + num_request_ranks;
    yac_proc_sphere_part_do_bnd_circle_search(
      grid_pair->proc_sphere_part, result_bnd_circle,
      curr_request_ranks, num_requests + i);
 
    // remove requests for local process
    int new_num_requests = 0;
    for (int j = 0; j < num_requests[i]; ++j) {
      if (curr_request_ranks[j] == comm_rank) continue;
      if (new_num_requests != j)
        curr_request_ranks[new_num_requests] = curr_request_ranks[j];
      ++new_num_requests;
    }
 
    num_request_ranks += (size_t)(num_requests[i] = new_num_requests);
  }
 
  int * sendcounts, * recvcounts, * sdispls, * rdispls;
  yac_get_comm_buffers(
    1, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
 
  for (size_t i = 0; i < num_request_ranks; ++i) sendcounts[request_ranks[i]]++;
 
  yac_generate_alltoallv_args(
    1, sendcounts, recvcounts, sdispls, rdispls, comm);
  int num_src_msg = 0, num_dst_msg = 0;
  for (int i = 0; i < comm_size; ++i) {
    num_src_msg += (recvcounts[i] > 0);
    num_dst_msg += (sendcounts[i] > 0);
  }
 
  size_t recv_count =
    (size_t)(rdispls[comm_size-1]) + (size_t)(recvcounts[comm_size-1]);
 
  coordinate_pointer coord_buffer =
    xmalloc((num_request_ranks + recv_count) * sizeof(*coord_buffer));
  coordinate_pointer send_request_coords = coord_buffer;
  coordinate_pointer recv_request_coords = coord_buffer + num_request_ranks;
 
  // pack bounding circles
  for (size_t i = 0, k = 0; i < count; ++i)
    for (size_t j = 0; j < num_requests[i]; ++j, ++k)
      memcpy(send_request_coords[sdispls[request_ranks[k]+1]++],
             search_coords[i], 3 * sizeof(double));
 
  MPI_Datatype coord_dt = yac_get_coordinate_mpi_datatype(comm);
 
  // exchange requests to other processes
  yac_alltoallv_p2p(
    send_request_coords, sendcounts, sdispls,
    recv_request_coords, recvcounts, rdispls,
    sizeof(*send_request_coords), coord_dt, comm);
 
  yac_mpi_call(MPI_Type_free(&coord_dt), comm);
 
  num_results = xrealloc(num_results, recv_count * sizeof(*num_results));
  memset(num_results, 0, recv_count * sizeof(*num_results));
 
  // do local search for requests
  yac_point_sphere_part_search_NNN(
    sphere_part, recv_count, recv_request_coords, n, &cos_angles,
    &cos_angles_array_size, &result_coords, &result_coords_array_size,
    &result_points, &result_points_array_size, num_results);
  yac_delete_point_sphere_part_search(sphere_part);
  free(coord_buffer);
 
  for (size_t i = 0; i < recv_count; ++i)
    if (max_num_results < num_results[i]) max_num_results = num_results[i];
  temp_results =
    xrealloc(temp_results, max_num_results * sizeof(*temp_results));
 
  // extract results (only the ones, that are local core points)
  struct nnn_search_result * request_results =
    xmalloc(n * (num_request_ranks + recv_count) * sizeof(*request_results));
  struct nnn_search_result * send_request_results = request_results;
  struct nnn_search_result * recv_request_results =
    request_results + n * recv_count;
 
  for (size_t i = 0, k = 0; i < recv_count; ++i) {
 
    struct nnn_search_result * curr_results = send_request_results + i * n;
    for (size_t j = 0; j < n; ++j, ++k) {
      size_t curr_local_id = result_points[k];
      curr_results[j].point_info.global_id = global_ids[curr_local_id];
      curr_results[j].point_info.data.rank = comm_rank;
      curr_results[j].point_info.data.orig_pos = (uint64_t)curr_local_id;
      curr_results[j].cos_angle = cos_angles[k];
    }
 
    for (size_t j = n; j < num_results[i]; ++j, ++k) {
      size_t curr_local_id = result_points[k];
      temp_results[j-n].point_info.global_id = global_ids[curr_local_id];
      temp_results[j-n].point_info.data.rank = comm_rank;
      temp_results[j-n].point_info.data.orig_pos = (uint64_t)curr_local_id;
      temp_results[j-n].cos_angle = cos_angles[k];
    }
    insert_nnn_result_points(curr_results, n, temp_results, num_results[i]-n);
  }
 
  free(cos_angles);
  free(result_coords);
  free(result_points);
  free(num_results);
 
  MPI_Datatype nnn_search_result_dt =
    yac_get_nnn_search_result_mpi_datatype(n, comm);
 
  // return results
  yac_alltoallv_p2p(
    send_request_results, recvcounts, rdispls,
    recv_request_results, sendcounts, sdispls,
    n * sizeof(*send_request_results), nnn_search_result_dt, comm);
 
  yac_mpi_call(MPI_Type_free(&nnn_search_result_dt), comm);
 
  // unpack results of remote search and integrate them into local results
  for (size_t i = 0, k = 0; i < count; ++i) {
    struct nnn_search_result * curr_results = results + n * i;
    for (size_t j = 0; j < num_requests[i]; ++k, ++j) {
      size_t pos = (size_t)(sdispls[request_ranks[k]]++);
      struct nnn_search_result * curr_request_results =
        recv_request_results + n * pos;
      insert_nnn_result_points(curr_results, n, curr_request_results, n);
    }
  }
  free(request_ranks);
  free(request_results);
  free(num_requests);
  yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
 
  size_t remote_count = 0;
  for (size_t i = 0; i < n * count; ++i)
    if (results[i].point_info.data.rank != comm_rank) ++remote_count;
 
  // gather remote points from search results
  struct single_remote_point * remote_points =
    xmalloc(remote_count * sizeof(*remote_points));
  for (size_t i = 0, k = 0; i < n * count; ++i)
    if (results[i].point_info.data.rank != comm_rank)
      remote_points[k++] = results[i].point_info;
 
  // convert all remote ids to local ones, extend local dist_grid data,
  // if necessary
  yac_dist_grid_single_remote_point_to_local(
    grid, remote_points, remote_count, field.location,
    local_ids + n * count - remote_count);
 
  free(remote_points);
 
  // gather final result
  for (size_t i = 0, offset = n * count - remote_count; i < n * count; ++i) {
    if (results[i].point_info.data.rank == comm_rank)
      local_ids[i] = (size_t)(results[i].point_info.data.orig_pos);
    else local_ids[i] = local_ids[offset++];
  }
  free(temp_results);
  free(results);
}
 
void yac_dist_grid_pair_do_bnd_circle_search(
  struct dist_grid_pair * grid_pair, int grid_id,
  const_bounding_circle_pointer bnd_circles, size_t count, size_t ** cells,
  size_t * num_results_per_bnd_circle, struct interp_field field) {
 
  MPI_Comm comm = grid_pair->comm;
  int comm_rank, comm_size;
  yac_mpi_call(MPI_Comm_rank(comm, &comm_rank), comm);
  yac_mpi_call(MPI_Comm_size(comm, &comm_size), comm);
 
  struct dist_grid * dist_grid =
    yac_dist_grid_pair_get_dist_grid(grid_pair, grid_id);
 
  int * num_ranks = xcalloc(count, sizeof(*num_ranks));
  int * rank_buffer = NULL;
  size_t rank_buffer_size = 0;
  size_t rank_buffer_array_size = 0;
 
  for (size_t i = 0; i < count; ++i) {
 
    ENSURE_ARRAY_SIZE(rank_buffer, rank_buffer_array_size,
                      rank_buffer_size + (size_t)comm_size);
 
    // finds all processes whose core area overlaps with the bounding circle
    // of the current cell
    // beware: even if the core are of a process does not overlap with the
    // bounding circle, it may have core cells that overlap nevertheless
    yac_proc_sphere_part_do_bnd_circle_search(
      grid_pair->proc_sphere_part, bnd_circles[i],
      rank_buffer + rank_buffer_size, num_ranks + i);
    rank_buffer_size += (size_t)(num_ranks[i]);
  }
 
  int * int_buffer =
    xmalloc(4 * (size_t)comm_size * sizeof(*int_buffer));
  int * result_sendcounts = int_buffer + 0 * comm_size;
  int * result_recvcounts = int_buffer + 1 * comm_size;
  int * result_sdispls =    int_buffer + 2 * comm_size;
  int * result_rdispls =    int_buffer + 3 * comm_size;
 
  int * sendcounts, * recvcounts, * sdispls, * rdispls;
  yac_get_comm_buffers(
    1, &sendcounts, &recvcounts, &sdispls, &rdispls, comm);
 
  for (size_t i = 0, offset = 0; i < count; ++i) {
    int curr_num_ranks = num_ranks[i];
    int * ranks = rank_buffer + offset;
    offset += (size_t)curr_num_ranks;
    for (int j = 0; j < curr_num_ranks; ++j) sendcounts[ranks[j]]++;
  }
 
  // local overlaps do not need to be send around
  size_t local_count = (size_t)(sendcounts[comm_rank]);
  sendcounts[comm_rank] = 0;
 
  yac_generate_alltoallv_args(
    1, sendcounts, recvcounts, sdispls, rdispls, comm);
 
  size_t send_count =
    (size_t)(sdispls[comm_size]) + (size_t)(sendcounts[comm_size-1]);
  size_t recv_count =
    (size_t)(rdispls[comm_size-1]) + (size_t)(recvcounts[comm_size-1]);
 
  struct bounding_circle * bnd_circle_buffer =
    xmalloc((send_count + recv_count + local_count) *
            sizeof(*bnd_circle_buffer));
  struct bounding_circle * send_buffer = bnd_circle_buffer;
  struct bounding_circle * recv_buffer = bnd_circle_buffer + send_count;
  struct bounding_circle * local_buffer =
    bnd_circle_buffer + send_count + recv_count;
 
  // pack bounding circles
  for (size_t i = 0, offset = 0, local_offset = 0; i < count; ++i) {
    int curr_num_ranks = num_ranks[i];
    int * ranks = rank_buffer + offset;
    offset += (size_t)curr_num_ranks;
    for (int j = 0; j < curr_num_ranks; ++j) {
      int rank = ranks[j];
      if (rank == comm_rank)
        local_buffer[local_offset++] = bnd_circles[i];
      else
        send_buffer[sdispls[rank + 1]++] = bnd_circles[i];
    }
  }
 
  MPI_Datatype bnd_circle_dt = yac_get_bounding_circle_mpi_datatype(comm);
  yac_mpi_call(MPI_Type_commit(&bnd_circle_dt), comm);
 
  // exchange bounding circles
  yac_alltoallv_p2p(
    send_buffer, sendcounts, sdispls, recv_buffer, recvcounts, rdispls,
    sizeof(*send_buffer), bnd_circle_dt, comm);
 
  yac_mpi_call(MPI_Type_free(&bnd_circle_dt), comm);
 
  struct bnd_sphere_part_search * cell_sphere_part =
    dist_grid_pair_get_cell_sphere_part(grid_pair, grid_id);
 
  size_t * local_cells = NULL;
  size_t * num_local_cells_per_bnd_circle =
    xmalloc((recv_count + local_count) *
            sizeof(*num_local_cells_per_bnd_circle));
 
  uint64_t * uint64_t_buffer =
    xmalloc((send_count + recv_count + local_count) *
            sizeof(*uint64_t_buffer));
  uint64_t * num_local_cells_per_bnd_circle_uint64_t = uint64_t_buffer;
  uint64_t * num_remote_cells_per_bnd_circle =
    uint64_t_buffer + recv_count + local_count;
 
  // search for received bounding circles in local data
  yac_bnd_sphere_part_search_do_bnd_circle_search(
    cell_sphere_part, recv_buffer, recv_count + local_count, &local_cells,
    num_local_cells_per_bnd_circle);
 
  int const * field_mask =
    (field.location == CELL)?
      yac_dist_grid_get_field_mask(dist_grid, field):NULL;
 
  struct bounding_circle * cell_bnd_circles = dist_grid->cell_bnd_circles;
 
  // check field mask and actual overlap of bounding circles
  for (size_t i = 0, offset = 0, new_offset = 0;
       i < recv_count + local_count; ++i) {
 
    struct bounding_circle * curr_bnd_circle = recv_buffer + i;
    size_t curr_num_results = num_local_cells_per_bnd_circle[i];
 
    // remove cells whose bounding circle do not overlap with the current one
    uint64_t new_num_results = 0;
    for (size_t j = 0; j < curr_num_results; ++j, ++offset) {
      size_t local_cell_id = local_cells[offset];
      if (!yac_extents_overlap(curr_bnd_circle,
                               cell_bnd_circles + local_cell_id)) continue;
      if ((field_mask == NULL) || (field_mask[local_cell_id])) {
        if (offset != new_offset) local_cells[new_offset] = local_cell_id;
        new_num_results++;
        new_offset++;
      }
    }
    num_local_cells_per_bnd_circle_uint64_t[i] = new_num_results;
  }
  free(num_local_cells_per_bnd_circle);
  free(bnd_circle_buffer);
 
  // exchange number of results per bounding circle
  yac_alltoallv_p2p(
    num_local_cells_per_bnd_circle_uint64_t, recvcounts, rdispls,
    num_remote_cells_per_bnd_circle, sendcounts, sdispls,
    sizeof(*num_local_cells_per_bnd_circle_uint64_t), MPI_UINT64_T, comm);
 
  for (int i = 0, saccu = 0, raccu = 0, soffset = 0, roffset = 0;
       i < comm_size; ++i) {
 
    result_sdispls[i] = saccu;
    result_rdispls[i] = raccu;
 
    int sendcount = recvcounts[i];
    int recvcount = sendcounts[i];
 
    result_sendcounts[i] = 0;
    result_recvcounts[i] = 0;
    for (int j = 0; j < sendcount; ++j, ++soffset)
      result_sendcounts[i] +=
        (int)(num_local_cells_per_bnd_circle_uint64_t[soffset]);
    for (int j = 0; j < recvcount; ++j, ++roffset)
      result_recvcounts[i] += (int)(num_remote_cells_per_bnd_circle[roffset]);
 
    saccu += result_sendcounts[i];
    raccu += result_recvcounts[i];
  }
 
  // count the number of results for bounding circles, which had a match with#
  // their original process
  size_t result_local_count = 0;
  for (size_t i = recv_count; i < recv_count + local_count; ++i)
    result_local_count += (size_t)(num_local_cells_per_bnd_circle_uint64_t[i]);
 
  size_t result_send_count = (size_t)(result_sdispls[comm_size-1]) +
                             (size_t)(result_sendcounts[comm_size-1]);
  size_t result_recv_count = (size_t)(result_rdispls[comm_size-1]) +
                             (size_t)(result_recvcounts[comm_size-1]);
 
  struct single_remote_point * single_remote_point_buffer =
    xmalloc((result_recv_count + result_send_count) *
            sizeof(*single_remote_point_buffer));
  struct single_remote_point * id_send_buffer = single_remote_point_buffer;
  struct single_remote_point * id_recv_buffer = single_remote_point_buffer +
                                                result_send_count;
 
  yac_int * cell_ids = dist_grid->cell_ids;
 
  for (size_t i = 0; i < result_send_count; ++i) {
    size_t local_cell_id = local_cells[i];
    id_send_buffer[i].global_id = cell_ids[local_cell_id];
    id_send_buffer[i].data.rank = comm_rank;
    id_send_buffer[i].data.orig_pos = local_cell_id;
  }
 
  MPI_Datatype single_remote_point_dt =
    yac_get_single_remote_point_mpi_datatype(comm);
 
  // redistribute results (global ids of found source cells)
  yac_alltoallv_p2p(
    id_send_buffer, result_sendcounts, result_sdispls,
    id_recv_buffer, result_recvcounts, result_rdispls,
    sizeof(*id_send_buffer), single_remote_point_dt, comm);
 
  yac_mpi_call(MPI_Type_free(&single_remote_point_dt), comm);
 
  size_t * new_local_cells =
    xmalloc((result_recv_count + result_local_count) *
            sizeof(*new_local_cells));
 
  memcpy(new_local_cells + result_recv_count,
         local_cells + result_send_count,
         result_local_count * sizeof(*new_local_cells));
  free(local_cells);
 
  // convert all remote ids to local ones, extend local dist_grid data,
  // if necessary
  yac_dist_grid_single_remote_point_to_local(
    dist_grid, id_recv_buffer, result_recv_count, CELL, new_local_cells);
 
  free(single_remote_point_buffer);
 
  size_t * reorder_idx =
    xmalloc((result_recv_count + result_local_count) * sizeof(*reorder_idx));
 
  memset(
    num_results_per_bnd_circle, 0, count * sizeof(*num_results_per_bnd_circle));
 
  for (size_t i = 0, offset = 0, reorder = 0, local_search_idx = recv_count,
              local_offset = result_recv_count; i < count; ++i) {
    int curr_num_ranks = num_ranks[i];
    int * ranks = rank_buffer + offset;
    offset += (size_t)curr_num_ranks;
    for (int j = 0; j < curr_num_ranks; ++j) {
      int rank = ranks[j];
      if (rank == comm_rank) {
        uint64_t curr_num_results =
          num_local_cells_per_bnd_circle_uint64_t[local_search_idx++];
        num_results_per_bnd_circle[i] += (size_t)curr_num_results;
        for (uint64_t k = 0; k < curr_num_results; ++k, ++reorder)
          reorder_idx[local_offset++] = reorder;
      } else {
        int rank_pos = sdispls[rank]++;
        uint64_t curr_num_results = num_remote_cells_per_bnd_circle[rank_pos];
        num_results_per_bnd_circle[i] += (size_t)curr_num_results;
        for (uint64_t k = 0; k < curr_num_results; ++k, ++reorder)
          reorder_idx[result_rdispls[rank]++] = reorder;
      }
    }
  }
  free(uint64_t_buffer);
  free(num_ranks);
  free(rank_buffer);
  free(int_buffer);
  yac_free_comm_buffers(sendcounts, recvcounts, sdispls, rdispls);
 
  yac_quicksort_index_size_t_size_t(
    reorder_idx, result_recv_count + result_local_count, new_local_cells);
  free(reorder_idx);
 
  // remove duplicated results
  for (size_t i = 0, offset = 0, new_offset = 0; i < count; ++i) {
 
    size_t * curr_local_cells = new_local_cells + offset;
    size_t curr_num_results_per_bnd_circle = num_results_per_bnd_circle[i];
    size_t new_num_results_per_bnd_circle = 0;
    size_t prev_cell = SIZE_MAX;
    offset += curr_num_results_per_bnd_circle;
 
    yac_quicksort_index_size_t_size_t(
      curr_local_cells, curr_num_results_per_bnd_circle, NULL);
 
    for (size_t j = 0; j < curr_num_results_per_bnd_circle; ++j) {
      size_t curr_cell = curr_local_cells[j];
      if (curr_cell != prev_cell) {
        new_local_cells[new_offset++] = (prev_cell = curr_cell);
        ++new_num_results_per_bnd_circle;
      }
    }
    num_results_per_bnd_circle[i] = new_num_results_per_bnd_circle;
  }
 
  *cells = new_local_cells;
}
 
void yac_dist_grid_pair_do_cell_search(
  struct dist_grid_pair * grid_pair, int search_grid_id, int result_grid_id,
  size_t * search_cells, size_t count, size_t ** result_cells,
  size_t * num_results_per_search_cell, struct