grid_tools.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2001 - 2018 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_grid_tools_h
#define dealii_grid_tools_h


#include <deal.II/base/bounding_box.h>
#include <deal.II/base/config.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/dofs/dof_handler.h>
#include <deal.II/fe/mapping.h>
#include <deal.II/fe/mapping_q1.h>
#include <deal.II/grid/manifold.h>
#include <deal.II/grid/tria.h>
#include <deal.II/grid/tria_accessor.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/hp/dof_handler.h>
#include <deal.II/lac/sparsity_tools.h>

#include <boost/optional.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/archive/binary_iarchive.hpp>
#include <boost/serialization/vector.hpp>
#include <boost/serialization/array.hpp>

#ifdef DEAL_II_WITH_ZLIB
#  include <boost/iostreams/stream.hpp>
#  include <boost/iostreams/filtering_stream.hpp>
#  include <boost/iostreams/device/back_inserter.hpp>
#  include <boost/iostreams/filter/gzip.hpp>
#endif


#include <bitset>
#include <list>
#include <set>

DEAL_II_NAMESPACE_OPEN

namespace parallel
{
  namespace distributed
  {
    template <int, int> class Triangulation;
  }
}

namespace hp
{
  template <int, int> class MappingCollection;
}

class SparsityPattern;

namespace internal
{
  template <int dim, int spacedim, class MeshType>
  class ActiveCellIterator
  {
  public:
#ifndef _MSC_VER
    typedef typename MeshType::active_cell_iterator type;
#else
    typedef TriaActiveIterator<::CellAccessor<dim, spacedim> > type;
#endif
  };

#ifdef _MSC_VER
  template <int dim, int spacedim>
  class ActiveCellIterator<dim, spacedim, ::DoFHandler<dim, spacedim> >
  {
  public:
    typedef TriaActiveIterator < ::DoFCellAccessor < ::DoFHandler<dim, spacedim >, false> > type;
  };

  template <int dim, int spacedim>
  class ActiveCellIterator<dim, spacedim, ::hp::DoFHandler<dim, spacedim> >
  {
  public:
    typedef TriaActiveIterator < ::DoFCellAccessor < ::hp::DoFHandler<dim, spacedim >, false> > type;
  };
#endif
}

namespace GridTools
{
  template <int dim, int spacedim>
  class Cache;


  template <int dim, int spacedim>
  double diameter (const Triangulation<dim, spacedim> &tria);

  template <int dim, int spacedim>
  double volume (const Triangulation<dim,spacedim> &tria,
                 const Mapping<dim,spacedim> &mapping = (StaticMappingQ1<dim,spacedim>::mapping));

  template <int dim, int spacedim>
  double
  minimal_cell_diameter (const Triangulation<dim, spacedim> &triangulation,
                         const Mapping<dim,spacedim> &mapping = (StaticMappingQ1<dim,spacedim>::mapping));

  template <int dim, int spacedim>
  double
  maximal_cell_diameter (const Triangulation<dim, spacedim> &triangulation,
                         const Mapping<dim,spacedim> &mapping = (StaticMappingQ1<dim,spacedim>::mapping));

  template <int dim>
  double cell_measure (const std::vector<Point<dim> > &all_vertices,
                       const unsigned int (&vertex_indices)[GeometryInfo<dim>::vertices_per_cell]);

  template <int dim, typename T>
  double cell_measure (const T &, ...);

  template <int dim, int spacedim>
  BoundingBox<spacedim> compute_bounding_box(const Triangulation<dim, spacedim> &triangulation);

  template <typename Iterator>
  Point<Iterator::AccessorType::space_dimension>
  project_to_object(const Iterator &object,
                    const Point<Iterator::AccessorType::space_dimension> &trial_point);


  template <int dim, int spacedim>
  void delete_unused_vertices (std::vector<Point<spacedim> >    &vertices,
                               std::vector<CellData<dim> > &cells,
                               SubCellData                 &subcelldata);

  template <int dim, int spacedim>
  void delete_duplicated_vertices (std::vector<Point<spacedim> >    &all_vertices,
                                   std::vector<CellData<dim> > &cells,
                                   SubCellData                 &subcelldata,
                                   std::vector<unsigned int>   &considered_vertices,
                                   const double                 tol=1e-12);


  template <int dim, typename Transformation, int spacedim>
  void transform (const Transformation        &transformation,
                  Triangulation<dim,spacedim> &triangulation);

  template <int dim, int spacedim>
  void shift (const Tensor<1,spacedim>    &shift_vector,
              Triangulation<dim,spacedim> &triangulation);


  void rotate (const double      angle,
               Triangulation<2> &triangulation);

  template <int dim>
  void
  rotate (const double          angle,
          const unsigned int    axis,
          Triangulation<dim,3> &triangulation);

  template <int dim>
  void laplace_transform (const std::map<unsigned int,Point<dim> > &new_points,
                          Triangulation<dim> &tria,
                          const Function<dim,double> *coefficient = nullptr,
                          const bool solve_for_absolute_positions = false);

  template <int dim, int spacedim>
  std::map<unsigned int,Point<spacedim> >
  get_all_vertices_at_boundary (const Triangulation<dim, spacedim> &tria);

  template <int dim, int spacedim>
  void scale (const double        scaling_factor,
              Triangulation<dim, spacedim> &triangulation);

  template <int dim, int spacedim>
  void distort_random (const double factor,
                       Triangulation<dim, spacedim> &triangulation,
                       const bool   keep_boundary=true);

  template <int dim, int spacedim>
  void
  remove_hanging_nodes( Triangulation<dim,spacedim> &tria,
                        const bool isotropic = false,
                        const unsigned int max_iterations = 100);

  template <int dim, int spacedim>
  void
  remove_anisotropy(  Triangulation<dim,spacedim> &tria,
                      const double max_ratio = 1.6180339887,
                      const unsigned int max_iterations = 5);

  template <int dim, int spacedim>
  void
  regularize_corner_cells(Triangulation<dim,spacedim> &tria,
                          const double limit_angle_fraction=.75);


  template <int dim, int spacedim>
#ifndef DOXYGEN
  std::tuple<
  std::vector<typename Triangulation<dim, spacedim>::active_cell_iterator >,
      std::vector< std::vector< Point<dim> > >,
      std::vector< std::vector<unsigned int> > >
#else
  return_type
#endif
      compute_point_locations(const Cache<dim,spacedim>                                         &cache,
                              const std::vector<Point<spacedim> >                               &points,
                              const typename Triangulation<dim, spacedim>::active_cell_iterator &cell_hint
                              = typename Triangulation<dim, spacedim>::active_cell_iterator());

  template <int dim, int spacedim>
#ifndef DOXYGEN
  std::tuple<
  std::vector< typename Triangulation<dim, spacedim>::active_cell_iterator >,
      std::vector< std::vector< Point<dim> > >,
      std::vector< std::vector< unsigned int > >,
      std::vector< std::vector< Point<spacedim> > >,
      std::vector< std::vector< unsigned int > >
      >
#else
  return_type
#endif
      distributed_compute_point_locations
      (const GridTools::Cache<dim,spacedim>                        &cache,
       const std::vector<Point<spacedim> >                         &local_points,
       const std::vector< std::vector< BoundingBox<spacedim> > >   &global_bboxes);

  template <int dim, int spacedim>
  std::map<unsigned int,Point<spacedim>> extract_used_vertices (
                                        const Triangulation<dim,spacedim> &container,
                                        const Mapping<dim,spacedim> &mapping = StaticMappingQ1<dim,spacedim>::mapping);

  template<int spacedim>
  unsigned int
  find_closest_vertex (const std::map<unsigned int,Point<spacedim>> &vertices,
                       const Point<spacedim>         &p);

  template <int dim, template <int, int> class MeshType, int spacedim>
  unsigned int
  find_closest_vertex (const MeshType<dim, spacedim> &mesh,
                       const Point<spacedim>         &p,
                       const std::vector<bool>       &marked_vertices = std::vector<bool>());

  template <int dim, template <int, int> class MeshType, int spacedim>
  unsigned int
  find_closest_vertex (const Mapping<dim,spacedim>   &mapping,
                       const MeshType<dim, spacedim> &mesh,
                       const Point<spacedim>         &p,
                       const std::vector<bool>       &marked_vertices = std::vector<bool>());


  template <int dim, template <int, int> class MeshType, int spacedim>
#ifndef _MSC_VER
  std::vector<typename MeshType<dim, spacedim>::active_cell_iterator>
#else
  std::vector<typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type>
#endif
  find_cells_adjacent_to_vertex (const MeshType<dim,spacedim> &container,
                                 const unsigned int            vertex_index);


  template <int dim, template <int,int> class MeshType, int spacedim>
#ifndef _MSC_VER
  typename MeshType<dim,spacedim>::active_cell_iterator
#else
  typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type
#endif
  find_active_cell_around_point (const MeshType<dim,spacedim> &mesh,
                                 const Point<spacedim>        &p,
                                 const std::vector<bool>      &marked_vertices = std::vector<bool>());

  template <int dim, template <int, int> class MeshType, int spacedim>
#ifndef _MSC_VER
  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> >
#else
  std::pair<typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type, Point<dim> >
#endif
  find_active_cell_around_point (const Mapping<dim,spacedim>  &mapping,
                                 const MeshType<dim,spacedim> &mesh,
                                 const Point<spacedim>        &p,
                                 const std::vector<bool>      &marked_vertices = std::vector<bool>());

  template <int dim, template <int, int> class MeshType, int spacedim>
#ifndef _MSC_VER
  std::pair<typename MeshType<dim, spacedim>::active_cell_iterator, Point<dim> >
#else
  std::pair<typename ::internal::ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim> >::type, Point<dim> >
#endif
  find_active_cell_around_point (const Mapping<dim,spacedim>                                                          &mapping,
                                 const MeshType<dim,spacedim>                                                         &mesh,
                                 const Point<spacedim>                                                                &p,
                                 const std::vector<std::set<typename MeshType<dim,spacedim>::active_cell_iterator > > &vertex_to_cell_map,
                                 const std::vector<std::vector<Tensor<1,spacedim> > >                                 &vertex_to_cell_centers,
                                 const typename MeshType<dim, spacedim>::active_cell_iterator                         &cell_hint=typename MeshType<dim, spacedim>::active_cell_iterator(),
                                 const std::vector<bool>                                                              &marked_vertices = std::vector<bool>());

  template <int dim, int spacedim>
  std::pair<typename hp::DoFHandler<dim, spacedim>::active_cell_iterator, Point<dim> >
  find_active_cell_around_point (const hp::MappingCollection<dim,spacedim> &mapping,
                                 const hp::DoFHandler<dim,spacedim>        &mesh,
                                 const Point<spacedim>                     &p);

  template <int dim, int spacedim>
  std::pair<typename Triangulation<dim, spacedim>::active_cell_iterator, Point<dim> >
  find_active_cell_around_point (const Cache<dim,spacedim>                                         &cache,
                                 const Point<spacedim>                                             &p,
                                 const typename Triangulation<dim, spacedim>::active_cell_iterator &cell_hint=typename Triangulation<dim, spacedim>::active_cell_iterator(),
                                 const std::vector<bool>                                           &marked_vertices = std::vector<bool>());

  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  get_active_child_cells (const typename MeshType::cell_iterator &cell);

  template <class MeshType>
  void
  get_active_neighbors (const typename MeshType::active_cell_iterator        &cell,
                        std::vector<typename MeshType::active_cell_iterator> &active_neighbors);

  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  compute_active_cell_halo_layer
  (const MeshType                                                                    &mesh,
   const std::function<bool (const typename MeshType::active_cell_iterator &)> &predicate);


  template <class MeshType>
  std::vector<typename MeshType::cell_iterator>
  compute_cell_halo_layer_on_level
  (const MeshType                                                       &mesh,
   const std::function<bool (const typename MeshType::cell_iterator &)> &predicate,
   const unsigned int                                                    level);


  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  compute_ghost_cell_halo_layer (const MeshType &mesh);

  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  compute_active_cell_layer_within_distance
  (const MeshType                                                                    &mesh,
   const std::function<bool (const typename MeshType::active_cell_iterator &)> &predicate,
   const double                                                                       layer_thickness);

  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  compute_ghost_cell_layer_within_distance ( const MeshType &mesh,
                                             const double layer_thickness);

  template <class MeshType>
  std::pair< Point<MeshType::space_dimension>, Point<MeshType::space_dimension> >
  compute_bounding_box
  ( const MeshType                                                                    &mesh,
    const std::function<bool (const typename MeshType::active_cell_iterator &)> &predicate );

  template < class MeshType >
  std::vector< BoundingBox< MeshType::space_dimension > >
  compute_mesh_predicate_bounding_box
  ( const MeshType                                                              &mesh,
    const std::function<bool (const typename MeshType::active_cell_iterator &)> &predicate,
    const unsigned int                                                           refinement_level = 0,
    const bool                                                                   allow_merge = false,
    const unsigned int                                                           max_boxes = numbers::invalid_unsigned_int);

  template <int spacedim>
#ifndef DOXYGEN
  std::tuple< std::vector< std::vector< unsigned int > >,
      std::map< unsigned int, unsigned int>,
      std::map< unsigned int, std::vector< unsigned int > > >
#else
  return_type
#endif
      guess_point_owner (const std::vector< std::vector< BoundingBox<spacedim> > >
                         &global_bboxes,
                         const std::vector< Point<spacedim> >    &points);


  template <int dim, int spacedim>
  std::vector<std::set<typename Triangulation<dim,spacedim>::active_cell_iterator> >
  vertex_to_cell_map(const Triangulation<dim,spacedim> &triangulation);

  template <int dim, int spacedim>
  std::vector<std::vector<Tensor<1,spacedim> > >
  vertex_to_cell_centers_directions(const Triangulation<dim,spacedim> &mesh,
                                    const std::vector<std::set<typename Triangulation<dim,spacedim>::active_cell_iterator> > &vertex_to_cells);


  template <int dim, int spacedim>
  unsigned int
  find_closest_vertex_of_cell(const typename Triangulation<dim,spacedim>::active_cell_iterator &cell,
                              const Point<spacedim> &position);

  template <int dim, int spacedim>
  std::map<unsigned int, types::global_vertex_index>
  compute_local_to_global_vertex_index_map(
    const parallel::distributed::Triangulation<dim,spacedim> &triangulation);

  template <int dim, int spacedim>
  std::pair<unsigned int, double>
  get_longest_direction(typename Triangulation<dim, spacedim>::active_cell_iterator cell);


  template <int dim, int spacedim>
  void
  get_face_connectivity_of_cells (const Triangulation<dim, spacedim> &triangulation,
                                  DynamicSparsityPattern             &connectivity);

  template <int dim, int spacedim>
  void
  get_vertex_connectivity_of_cells (const Triangulation<dim, spacedim> &triangulation,
                                    DynamicSparsityPattern             &connectivity);

  template <int dim, int spacedim>
  void
  get_vertex_connectivity_of_cells_on_level (const Triangulation<dim, spacedim> &triangulation,
                                             const unsigned int                 level,
                                             DynamicSparsityPattern             &connectivity);

  template <int dim, int spacedim>
  void
  partition_triangulation (const unsigned int  n_partitions,
                           Triangulation<dim, spacedim> &triangulation,
                           const SparsityTools::Partitioner partitioner = SparsityTools::Partitioner::metis
                          );

  template <int dim, int spacedim>
  void
  partition_triangulation (const unsigned int                n_partitions,
                           const std::vector<unsigned int>  &cell_weights,
                           Triangulation<dim, spacedim>     &triangulation,
                           const SparsityTools::Partitioner  partitioner = SparsityTools::Partitioner::metis
                          );

  template <int dim, int spacedim>
  void
  partition_triangulation (const unsigned int     n_partitions,
                           const SparsityPattern &cell_connection_graph,
                           Triangulation<dim,spacedim>    &triangulation,
                           const SparsityTools::Partitioner partitioner = SparsityTools::Partitioner::metis
                          );

  template <int dim, int spacedim>
  void
  partition_triangulation (const unsigned int                n_partitions,
                           const std::vector<unsigned int>  &cell_weights,
                           const SparsityPattern            &cell_connection_graph,
                           Triangulation<dim,spacedim>      &triangulation,
                           const SparsityTools::Partitioner  partitioner = SparsityTools::Partitioner::metis
                          );

  template <int dim, int spacedim>
  void
  partition_triangulation_zorder (const unsigned int          n_partitions,
                                  Triangulation<dim,spacedim> &triangulation);

  template <int dim, int spacedim>
  void
  partition_multigrid_levels (Triangulation<dim,spacedim> &triangulation);

  template <int dim, int spacedim>
  void
  get_subdomain_association (const Triangulation<dim, spacedim>  &triangulation,
                             std::vector<types::subdomain_id> &subdomain);

  template <int dim, int spacedim>
  unsigned int
  count_cells_with_subdomain_association (const Triangulation<dim, spacedim> &triangulation,
                                          const types::subdomain_id         subdomain);


  template <int dim, int spacedim>
  std::vector<bool>
  get_locally_owned_vertices (const Triangulation<dim,spacedim> &triangulation);


  template <typename MeshType>
  std::list<std::pair<typename MeshType::cell_iterator,
      typename MeshType::cell_iterator> >
      get_finest_common_cells (const MeshType &mesh_1,
                               const MeshType &mesh_2);

  template <int dim, int spacedim>
  bool
  have_same_coarse_mesh (const Triangulation<dim, spacedim> &mesh_1,
                         const Triangulation<dim, spacedim> &mesh_2);

  template <typename MeshType>
  bool
  have_same_coarse_mesh (const MeshType &mesh_1,
                         const MeshType &mesh_2);


  template <int dim, int spacedim>
  typename Triangulation<dim,spacedim>::DistortedCellList
  fix_up_distorted_child_cells (const typename Triangulation<dim,spacedim>::DistortedCellList &distorted_cells,
                                Triangulation<dim,spacedim> &triangulation);






  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  get_patch_around_cell(const typename MeshType::active_cell_iterator &cell);


  template <class Container>
  std::vector<typename Container::cell_iterator>
  get_cells_at_coarsest_common_level(const std::vector<typename Container::active_cell_iterator> &patch_cells);

  template <class Container>
  void
  build_triangulation_from_patch (
    const std::vector<typename Container::active_cell_iterator>  &patch,
    Triangulation<Container::dimension,Container::space_dimension> &local_triangulation,
    std::map<typename Triangulation<Container::dimension,Container::space_dimension>::active_cell_iterator,
    typename Container::active_cell_iterator> &patch_to_global_tria_map);

  template <class DoFHandlerType>
  std::map< types::global_dof_index,std::vector<typename DoFHandlerType::active_cell_iterator> >
  get_dof_to_support_patch_map(DoFHandlerType &dof_handler);



  template <typename CellIterator>
  struct PeriodicFacePair
  {
    CellIterator cell[2];

    unsigned int face_idx[2];

    std::bitset<3> orientation;

    FullMatrix<double> matrix;
  };


  template <typename FaceIterator>
  bool
  orthogonal_equality (std::bitset<3>     &orientation,
                       const FaceIterator &face1,
                       const FaceIterator &face2,
                       const int          direction,
                       const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset
                       = Tensor<1,FaceIterator::AccessorType::space_dimension>(),
                       const FullMatrix<double> &matrix = FullMatrix<double>());


  template <typename FaceIterator>
  bool
  orthogonal_equality (const FaceIterator &face1,
                       const FaceIterator &face2,
                       const int          direction,
                       const Tensor<1,FaceIterator::AccessorType::space_dimension> &offset
                       = Tensor<1,FaceIterator::AccessorType::space_dimension>(),
                       const FullMatrix<double> &matrix = FullMatrix<double>());


  template <typename MeshType>
  void
  collect_periodic_faces
  (const MeshType                            &mesh,
   const types::boundary_id                   b_id1,
   const types::boundary_id                   b_id2,
   const int                                  direction,
   std::vector<PeriodicFacePair<typename MeshType::cell_iterator> > &matched_pairs,
   const Tensor<1,MeshType::space_dimension> &offset = ::Tensor<1,MeshType::space_dimension>(),
   const FullMatrix<double>                  &matrix = FullMatrix<double>());


  template <typename MeshType>
  void
  collect_periodic_faces
  (const MeshType                                    &mesh,
   const types::boundary_id                           b_id,
   const int                                          direction,
   std::vector<PeriodicFacePair<typename MeshType::cell_iterator> > &matched_pairs,
   const ::Tensor<1,MeshType::space_dimension> &offset = ::Tensor<1,MeshType::space_dimension>(),
   const FullMatrix<double>                          &matrix = FullMatrix<double>());


  template <int dim, int spacedim>
  void copy_boundary_to_manifold_id(Triangulation<dim, spacedim> &tria,
                                    const bool reset_boundary_ids=false);

  template <int dim, int spacedim>
  void map_boundary_to_manifold_ids(const std::vector<types::boundary_id> &src_boundary_ids,
                                    const std::vector<types::manifold_id> &dst_manifold_ids,
                                    Triangulation<dim, spacedim> &tria,
                                    const std::vector<types::boundary_id> &reset_boundary_ids = {});

  template <int dim, int spacedim>
  void copy_material_to_manifold_id(Triangulation<dim, spacedim> &tria,
                                    const bool compute_face_ids=false);


  template <typename DataType, typename MeshType>
  void
  exchange_cell_data_to_ghosts (const MeshType &mesh,
                                const std::function<boost::optional<DataType> (const typename MeshType::active_cell_iterator &)> &pack,
                                const std::function<void (const typename MeshType::active_cell_iterator &, const DataType &)> &unpack);

  /* Exchange with all processors of the MPI communicator @p mpi_communicator the vector of bounding
   * boxes @p local_bboxes.
   *
   * This function is meant to exchange bounding boxes describing the locally owned
   * cells in a distributed triangulation obtained with the function
   * GridTools::compute_mesh_predicate_bounding_box .
   *
   * The output vector's size is the number of processes of the MPI communicator:
   * its i-th entry contains the vector @p local_bboxes of the i-th process.
   */
  template<int spacedim>
  std::vector< std::vector< BoundingBox<spacedim> > >
  exchange_local_bounding_boxes(const std::vector< BoundingBox<spacedim> > &local_bboxes,
                                MPI_Comm                                    mpi_communicator);

  template <int dim, typename T>
  struct CellDataTransferBuffer
  {
    std::vector<CellId> cell_ids;

    std::vector<T> data;

    template <class Archive>
    void save (Archive &ar,
               const unsigned int version) const;

    template <class Archive>
    void load (Archive &ar,
               const unsigned int version);

    BOOST_SERIALIZATION_SPLIT_MEMBER()
  };


  DeclException1 (ExcInvalidNumberOfPartitions,
                  int,
                  << "The number of partitions you gave is " << arg1
                  << ", but must be greater than zero.");
  DeclException1 (ExcNonExistentSubdomain,
                  int,
                  << "The subdomain id " << arg1
                  << " has no cells associated with it.");
  DeclException0 (ExcTriangulationHasBeenRefined);

  DeclException1 (ExcScalingFactorNotPositive,
                  double,
                  << "The scaling factor must be positive, but it is " << arg1 << ".");
  template <int N>
  DeclException1 (ExcPointNotFoundInCoarseGrid,
                  Point<N>,
                  << "The point <" << arg1
                  << "> could not be found inside any of the "
                  << "coarse grid cells.");
  template <int N>
  DeclException1 (ExcPointNotFound,
                  Point<N>,
                  << "The point <" << arg1
                  << "> could not be found inside any of the "
                  << "subcells of a coarse grid cell.");

  DeclException1 (ExcVertexNotUsed,
                  unsigned int,
                  << "The given vertex with index " << arg1
                  << " is not used in the given triangulation.");


} /*namespace GridTools*/



/* ----------------- Template function --------------- */

#ifndef DOXYGEN

namespace GridTools
{
  template <int dim, typename T>
  double cell_measure (const T &, ...)
  {
    Assert(false, ExcNotImplemented());
    return std::numeric_limits<double>::quiet_NaN();
  }

  template <int dim, typename Predicate, int spacedim>
  void transform (const Predicate    &predicate,
                  Triangulation<dim, spacedim> &triangulation)
  {
    std::vector<bool> treated_vertices (triangulation.n_vertices(),
                                        false);

    // loop over all active cells, and
    // transform those vertices that
    // have not yet been touched. note
    // that we get to all vertices in
    // the triangulation by only
    // visiting the active cells.
    typename Triangulation<dim, spacedim>::active_cell_iterator
    cell = triangulation.begin_active (),
    endc = triangulation.end ();
    for (; cell!=endc; ++cell)
      for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
        if (treated_vertices[cell->vertex_index(v)] == false)
          {
            // transform this vertex
            cell->vertex(v) = predicate(cell->vertex(v));
            // and mark it as treated
            treated_vertices[cell->vertex_index(v)] = true;
          };


    // now fix any vertices on hanging nodes so that we don't create any holes
    if (dim==2)
      {
        typename Triangulation<dim,spacedim>::active_cell_iterator
        cell = triangulation.begin_active(),
        endc = triangulation.end();
        for (; cell!=endc; ++cell)
          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
            if (cell->face(face)->has_children() &&
                !cell->face(face)->at_boundary())
              {
                // this line has children
                cell->face(face)->child(0)->vertex(1)
                  = (cell->face(face)->vertex(0) +
                     cell->face(face)->vertex(1)) / 2;
              }
      }
    else if (dim==3)
      {
        typename Triangulation<dim,spacedim>::active_cell_iterator
        cell = triangulation.begin_active(),
        endc = triangulation.end();
        for (; cell!=endc; ++cell)
          for (unsigned int face=0; face<GeometryInfo<dim>::faces_per_cell; ++face)
            if (cell->face(face)->has_children() &&
                !cell->face(face)->at_boundary())
              {
                // this face has hanging nodes
                cell->face(face)->child(0)->vertex(1)
                  = (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) / 2.0;
                cell->face(face)->child(0)->vertex(2)
                  = (cell->face(face)->vertex(0) + cell->face(face)->vertex(2)) / 2.0;
                cell->face(face)->child(1)->vertex(3)
                  = (cell->face(face)->vertex(1) + cell->face(face)->vertex(3)) / 2.0;
                cell->face(face)->child(2)->vertex(3)
                  = (cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) / 2.0;

                // center of the face
                cell->face(face)->child(0)->vertex(3)
                  = (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)
                     + cell->face(face)->vertex(2) + cell->face(face)->vertex(3)) / 4.0;
              }
      }

    // Make sure FEValues notices that the mesh has changed
    triangulation.signals.mesh_movement();
  }



  template <class MeshType>
  std::vector<typename MeshType::active_cell_iterator>
  get_active_child_cells (const typename MeshType::cell_iterator &cell)
  {
    std::vector<typename MeshType::active_cell_iterator> child_cells;

    if (cell->has_children())
      {
        for (unsigned int child=0;
             child<cell->n_children(); ++child)
          if (cell->child (child)->has_children())
            {
              const std::vector<typename MeshType::active_cell_iterator>
              children = get_active_child_cells<MeshType> (cell->child(child));
              child_cells.insert (child_cells.end(),
                                  children.begin(), children.end());
            }
          else
            child_cells.push_back (cell->child(child));
      }

    return child_cells;
  }



  template <class MeshType>
  void
  get_active_neighbors(const typename MeshType::active_cell_iterator        &cell,
                       std::vector<typename MeshType::active_cell_iterator> &active_neighbors)
  {
    active_neighbors.clear ();
    for (unsigned int n=0; n<GeometryInfo<MeshType::dimension>::faces_per_cell; ++n)
      if (! cell->at_boundary(n))
        {
          if (MeshType::dimension == 1)
            {
              // check children of neighbor. note
              // that in 1d children of the neighbor
              // may be further refined. In 1d the
              // case is simple since we know what
              // children bound to the present cell
              typename MeshType::cell_iterator
              neighbor_child = cell->neighbor(n);
              if (!neighbor_child->active())
                {
                  while (neighbor_child->has_children())
                    neighbor_child = neighbor_child->child (n==0 ? 1 : 0);

                  Assert (neighbor_child->neighbor(n==0 ? 1 : 0)==cell,
                          ExcInternalError());
                }
              active_neighbors.push_back (neighbor_child);
            }
          else
            {
              if (cell->face(n)->has_children())
                // this neighbor has children. find
                // out which border to the present
                // cell
                for (unsigned int c=0; c<cell->face(n)->number_of_children(); ++c)
                  active_neighbors.push_back (cell->neighbor_child_on_subface(n,c));
              else
                {
                  // the neighbor must be active
                  // himself
                  Assert(cell->neighbor(n)->active(), ExcInternalError());
                  active_neighbors.push_back(cell->neighbor(n));
                }
            }
        }
  }



  namespace internal
  {
    namespace ProjectToObject
    {
      struct CrossDerivative
      {
        const unsigned int direction_0;
        const unsigned int direction_1;

        CrossDerivative(const unsigned int d0, const unsigned int d1);
      };

      inline
      CrossDerivative::CrossDerivative(const unsigned int d0, const unsigned int d1)
        :
        direction_0 (d0),
        direction_1 (d1)
      {}



      template <typename F>
      inline
      auto
      centered_first_difference(const double  center,
                                const double  step,
                                const F      &f)
      -> decltype(f(center) - f(center))
      {
        return (f(center + step) - f(center - step))/(2.0*step);
      }



      template <typename F>
      inline
      auto
      centered_second_difference(const double  center,
                                 const double  step,
                                 const F      &f)
      -> decltype(f(center) - f(center))
      {
        return (f(center + step) - 2.0*f(center) + f(center - step))/(step*step);
      }



      template <int structdim, typename F>
      inline
      auto
      cross_stencil
      (const CrossDerivative                                        cross_derivative,
       const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &center,
       const double                                                 step,
       const F                                                     &f)
      -> decltype(f(center) - f(center))
      {
        Tensor<1, GeometryInfo<structdim>::vertices_per_cell> simplex_vector;
        simplex_vector[cross_derivative.direction_0] = 0.5*step;
        simplex_vector[cross_derivative.direction_1] = -0.5*step;
        return (- 4.0     *f(center)
                - 1.0     *f(center + simplex_vector)
                - 1.0/3.0 *f(center - simplex_vector)
                + 16.0/3.0*f(center + 0.5*simplex_vector)
               )/step;
      }



      template <int spacedim, int structdim, typename F>
      inline
      double
      gradient_entry
      (const unsigned int                                           row_n,
       const unsigned int                                           dependent_direction,
       const Point<spacedim>                                       &p0,
       const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &center,
       const double                                                 step,
       const F                                                     &f)
      {
        Assert(row_n < GeometryInfo<structdim>::vertices_per_cell &&
               dependent_direction < GeometryInfo<structdim>::vertices_per_cell,
               ExcMessage("This function assumes that the last weight is a "
                          "dependent variable (and hence we cannot take its "
                          "derivative directly)."));
        Assert(row_n != dependent_direction,
               ExcMessage("We cannot differentiate with respect to the variable "
                          "that is assumed to be dependent."));

        const Point<spacedim> manifold_point = f(center);
        const Tensor<1, spacedim> stencil_value = cross_stencil<structdim>
                                                  ({row_n, dependent_direction},
                                                   center,
                                                   step,
                                                   f);
        double entry = 0.0;
        for (unsigned int dim_n = 0; dim_n < spacedim; ++dim_n)
          entry += -2.0*(p0[dim_n] - manifold_point[dim_n])*stencil_value[dim_n];
        return entry;
      }

      template <typename Iterator, int spacedim, int structdim>
      Point<spacedim>
      project_to_d_linear_object (const Iterator        &object,
                                  const Point<spacedim> &trial_point)
      {
        // let's look at this for simplicity for a quad (structdim==2) in a space with
        // spacedim>2 (notate trial_point by y): all points on the surface are
        // given by
        //   x(\xi) = sum_i v_i phi_x(\xi)
        // where v_i are the vertices of the quad, and \xi=(\xi_1,\xi_2) are the
        // reference coordinates of the quad. so what we are trying to do is find
        // a point x on the surface that is closest to the point y. there are
        // different ways to solve this problem, but in the end it's a nonlinear
        // problem and we have to find reference coordinates \xi so that J(\xi) =
        // 1/2 || x(\xi)-y ||^2 is minimal. x(\xi) is a function that is
        // structdim-linear in \xi, so J(\xi) is a polynomial of degree 2*structdim that we'd
        // like to minimize. unless structdim==1, we'll have to use a Newton method to
        // find the answer. This leads to the following formulation of Newton
        // steps:
        //
        // Given \xi_k, find \delta\xi_k so that
        //   H_k \delta\xi_k = - F_k
        // where H_k is an approximation to the second derivatives of J at \xi_k,
        // and F_k is the first derivative of J.  We'll iterate this a number of
        // times until the right hand side is small enough. As a stopping
        // criterion, we terminate if ||\delta\xi||<eps.
        //
        // As for the Hessian, the best choice would be
        //   H_k = J''(\xi_k)
        // but we'll opt for the simpler Gauss-Newton form
        //   H_k = A^T A
        // i.e.
        //   (H_k)_{nm} = \sum_{i,j} v_i*v_j *
        //                   \partial_n phi_i *
        //                   \partial_m phi_j
        // we start at xi=(0.5, 0.5).
        Point<structdim> xi;
        for (unsigned int d=0; d<structdim; ++d)
          xi[d] = 0.5;

        Point<spacedim> x_k;
        for (unsigned int i=0; i<GeometryInfo<structdim>::vertices_per_cell; ++i)
          x_k += object->vertex(i) *
                 GeometryInfo<structdim>::d_linear_shape_function (xi, i);

        do
          {
            Tensor<1,structdim> F_k;
            for (unsigned int i=0; i<GeometryInfo<structdim>::vertices_per_cell; ++i)
              F_k += (x_k-trial_point)*object->vertex(i) *
                     GeometryInfo<structdim>::d_linear_shape_function_gradient (xi, i);

            Tensor<2,structdim> H_k;
            for (unsigned int i=0; i<GeometryInfo<structdim>::vertices_per_cell; ++i)
              for (unsigned int j=0; j<GeometryInfo<structdim>::vertices_per_cell; ++j)
                {
                  Tensor<2, structdim> tmp = outer_product(
                                               GeometryInfo<structdim>::d_linear_shape_function_gradient(xi, i),
                                               GeometryInfo<structdim>::d_linear_shape_function_gradient(xi, j));
                  H_k += (object->vertex(i) * object->vertex(j)) * tmp;
                }

            const Tensor<1,structdim> delta_xi = - invert(H_k) * F_k;
            xi += delta_xi;

            x_k = Point<spacedim>();
            for (unsigned int i=0; i<GeometryInfo<structdim>::vertices_per_cell; ++i)
              x_k += object->vertex(i) *
                     GeometryInfo<structdim>::d_linear_shape_function (xi, i);

            if (delta_xi.norm() < 1e-7)
              break;
          }
        while (true);

        return x_k;
      }
    }
  }



  namespace internal
  {
    // We hit an internal compiler error in ICC 15 if we define this as a lambda
    // inside the project_to_object function below.
    template <int structdim>
    inline
    bool weights_are_ok (const Tensor<1, GeometryInfo<structdim>::vertices_per_cell> &v)
    {
      // clang has trouble figuring out structdim here, so define it
      // again:
      static const std::size_t n_vertices_per_cell
        = Tensor<1, GeometryInfo<structdim>::vertices_per_cell>::n_independent_components;
      std::array<double, n_vertices_per_cell> copied_weights;
      for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
        {
          copied_weights[i] = v[i];
          if (v[i] < 0.0 || v[i] > 1.0)
            return false;
        }

      // check the sum: try to avoid some roundoff errors by summing in order
      std::sort(copied_weights.begin(), copied_weights.end());
      const double sum = std::accumulate(copied_weights.begin(), copied_weights.end(), 0.0);
      return std::abs(sum - 1.0) < 1e-10; // same tolerance used in manifold.cc
    }
  }

  template <typename Iterator>
  Point<Iterator::AccessorType::space_dimension>
  project_to_object(const Iterator &object,
                    const Point<Iterator::AccessorType::space_dimension> &trial_point)
  {
    const int spacedim = Iterator::AccessorType::space_dimension;
    const int structdim = Iterator::AccessorType::structure_dimension;

    Point<spacedim> projected_point = trial_point;

    if (structdim >= spacedim)
      return projected_point;
    else if (structdim == 1 || structdim == 2)
      {
        using namespace internal::ProjectToObject;
        // Try to use the special flat algorithm for quads (this is better
        // than the general algorithm in 3D). This does not take into account
        // whether projected_point is outside the quad, but we optimize along
        // lines below anyway:
        const int dim = Iterator::AccessorType::dimension;
        const Manifold<dim, spacedim> &manifold = object->get_manifold();
        if (structdim == 2 &&
            dynamic_cast<const FlatManifold<dim,spacedim> *>(&manifold)
            != nullptr)
          {
            projected_point = project_to_d_linear_object<Iterator, spacedim, structdim>(object, trial_point);
          }
        else
          {
            // We want to find a point on the convex hull (defined by the
            // vertices of the object and the manifold description) that is
            // relatively close to the trial point. This has a few issues:
            //
            // 1. For a general convex hull we are not guaranteed that a unique
            //    minimum exists.
            // 2. The independent variables in the optimization process are the
            //    weights given to Manifold::get_new_point, which must sum to 1,
            //    so we cannot use standard finite differences to approximate a
            //    gradient.
            //
            // There is not much we can do about 1., but for 2. we can derive
            // finite difference stencils that work on a structdim-dimensional
            // simplex and rewrite the optimization problem to use those
            // instead. Consider the structdim 2 case and let
            //
            // F(c0, c1, c2, c3) = Manifold::get_new_point(vertices, {c0, c1, c2, c3})
            //
            // where {c0, c1, c2, c3} are the weights for the four vertices on
            // the quadrilateral. We seek to minimize the Euclidean distance
            // between F(...) and trial_point. We can solve for c3 in terms of
            // the other weights and get, for one coordinate direction
            //
            // d/dc0 ((x0 - F(c0, c1, c2, 1 - c0 - c1 - c2))^2)
            //      = -2(x0 - F(...)) (d/dc0 F(...) - d/dc3 F(...))
            //
            // where we substitute back in for c3 after taking the
            // derivative. We can compute a stencil for the cross derivative
            // d/dc0 - d/dc3: this is exactly what cross_stencil approximates
            // (and gradient_entry computes the sum over the independent
            // variables). Below, we somewhat arbitrarily pick the last
            // component as the dependent one.
            //
            // Since we can now calculate derivatives of the objective
            // function we can use gradient descent to minimize it.
            //
            // Of course, this is much simpler in the structdim = 1 case (we
            // could rewrite the projection as a 1D optimization problem), but
            // to reduce the potential for bugs we use the same code in both
            // cases.
            const double step_size = object->diameter()/64.0;

            constexpr unsigned int n_vertices_per_cell = GeometryInfo<structdim>::vertices_per_cell;

            std::array<Point<spacedim>, n_vertices_per_cell> vertices;
            for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
                 ++vertex_n)
              vertices[vertex_n] = object->vertex(vertex_n);

            auto get_point_from_weights =
              [&](const Tensor<1, n_vertices_per_cell> &weights)
              -> Point<spacedim>
            {
              return object->get_manifold().get_new_point
              (make_array_view(vertices.begin(), vertices.end()),
              make_array_view(&weights[0],
              &weights[n_vertices_per_cell - 1] + 1));
            };

            // pick the initial weights as (normalized) inverse distances from
            // the trial point:
            Tensor<1, n_vertices_per_cell> guess_weights;
            double guess_weights_sum = 0.0;
            for (unsigned int vertex_n = 0; vertex_n < n_vertices_per_cell;
                 ++vertex_n)
              {
                const double distance = vertices[vertex_n].distance(trial_point);
                if (distance == 0.0)
                  {
                    guess_weights = 0.0;
                    guess_weights[vertex_n] = 1.0;
                    guess_weights_sum = 1.0;
                    break;
                  }
                else
                  {
                    guess_weights[vertex_n] = 1.0/distance;
                    guess_weights_sum += guess_weights[vertex_n];
                  }
              }
            guess_weights /= guess_weights_sum;
            Assert(internal::weights_are_ok<structdim>(guess_weights), ExcInternalError());

            // The optimization algorithm consists of two parts:
            //
            // 1. An outer loop where we apply the gradient descent algorithm.
            // 2. An inner loop where we do a line search to find the optimal
            //    length of the step one should take in the gradient direction.
            //
            for (unsigned int outer_n = 0; outer_n < 40; ++outer_n)
              {
                const unsigned int dependent_direction = n_vertices_per_cell - 1;
                Tensor<1, n_vertices_per_cell> current_gradient;
                for (unsigned int row_n = 0;
                     row_n < n_vertices_per_cell;
                     ++row_n)
                  {
                    if (row_n != dependent_direction)
                      {
                        current_gradient[row_n] = gradient_entry<spacedim, structdim>
                                                  (row_n,
                                                   dependent_direction,
                                                   trial_point,
                                                   guess_weights,
                                                   step_size,
                                                   get_point_from_weights);

                        current_gradient[dependent_direction] -= current_gradient[row_n];
                      }
                  }

                // We need to travel in the -gradient direction, as noted
                // above, but we may not want to take a full step in that
                // direction; instead, guess that we will go -0.5*gradient and
                // do quasi-Newton iteration to pick the best multiplier. The
                // goal is to find a scalar alpha such that
                //
                // F(x - alpha g)
                //
                // is minimized, where g is the gradient and F is the
                // objective function. To find the optimal value we find roots
                // of the derivative of the objective function with respect to
                // alpha by Newton iteration, where we approximate the first
                // and second derivatives of F(x - alpha g) with centered
                // finite differences.
                double gradient_weight = -0.5;
                auto gradient_weight_objective_function = [&](const double gradient_weight_guess)
                                                          -> double
                {
                  return (trial_point -
                  get_point_from_weights(guess_weights +
                  gradient_weight_guess*current_gradient)).norm_square();
                };

                for (unsigned int inner_n = 0; inner_n < 10; ++inner_n)
                  {
                    const double update_numerator = centered_first_difference
                                                    (gradient_weight, step_size, gradient_weight_objective_function);
                    const double update_denominator = centered_second_difference
                                                      (gradient_weight, step_size, gradient_weight_objective_function);

                    // avoid division by zero. Note that we limit the gradient weight below
                    if (std::abs(update_denominator) == 0.0)
                      break;
                    gradient_weight = gradient_weight - update_numerator/update_denominator;

                    // Put a fairly lenient bound on the largest possible
                    // gradient (things tend to be locally flat, so the gradient
                    // itself is usually small)
                    if (std::abs(gradient_weight) > 10)
                      {
                        gradient_weight = -10.0;
                        break;
                      }
                  }

                // It only makes sense to take convex combinations with weights
                // between zero and one. If the update takes us outside of this
                // region then rescale the update to stay within the region and
                // try again
                Tensor<1, n_vertices_per_cell> tentative_weights =
                  guess_weights + gradient_weight*current_gradient;

                double new_gradient_weight = gradient_weight;
                for (unsigned int iteration_count = 0; iteration_count < 40; ++iteration_count)
                  {
                    if (internal::weights_are_ok<structdim>(tentative_weights))
                      break;

                    for (unsigned int i = 0; i < n_vertices_per_cell; ++i)
                      {
                        if (tentative_weights[i] < 0.0)
                          {
                            tentative_weights -= (tentative_weights[i]/current_gradient[i])
                                                 *current_gradient;
                          }
                        if (tentative_weights[i] < 0.0 || 1.0 < tentative_weights[i])
                          {
                            new_gradient_weight /= 2.0;
                            tentative_weights = guess_weights + new_gradient_weight*current_gradient;
                          }
                      }
                  }

                // the update might still send us outside the valid region, so
                // check again and quit if the update is still not valid
                if (!internal::weights_are_ok<structdim>(tentative_weights))
                  break;

                // if we cannot get closer by traveling in the gradient direction then quit
                if (get_point_from_weights(tentative_weights).distance(trial_point) <
                    get_point_from_weights(guess_weights).distance(trial_point))
                  guess_weights = tentative_weights;
                else
                  break;
                Assert(internal::weights_are_ok<structdim>(guess_weights), ExcInternalError());
              }
            Assert(internal::weights_are_ok<structdim>(guess_weights), ExcInternalError());
            projected_point =  get_point_from_weights(guess_weights);
          }

        // if structdim == 2 and the optimal point is not on the interior then
        // we may be able to get a more accurate result by projecting onto the
        // lines.
        if (structdim == 2)
          {
            std::array<Point<spacedim>, GeometryInfo<structdim>::lines_per_cell>
            line_projections;
            for (unsigned int line_n = 0; line_n < GeometryInfo<structdim>::lines_per_cell;
                 ++line_n)
              {
                line_projections[line_n] = project_to_object(object->line(line_n),
                                                             trial_point);
              }
            std::sort(line_projections.begin(), line_projections.end(),
                      [&](const Point<spacedim> &a, const Point<spacedim> &b)
            {
              return a.distance(trial_point) < b.distance(trial_point);
            });
            if (line_projections[0].distance(trial_point)
                < projected_point.distance(trial_point))
              projected_point = line_projections[0];
          }
      }
    else
      {
        Assert(false, ExcNotImplemented());
        return projected_point;
      }

    return projected_point;
  }



  template <int dim, typename T>
  template <class Archive>
  void
  CellDataTransferBuffer<dim,T>::save (Archive &ar,
                                       const unsigned int /*version*/) const
  {
    Assert(cell_ids.size() == data.size(),
           ExcDimensionMismatch(cell_ids.size(), data.size()));
    // archive the cellids in an efficient binary format
    const size_t n_cells = cell_ids.size();
    ar &n_cells;
    for (auto &it : cell_ids)
      {
        CellId::binary_type binary_cell_id = it.template to_binary<dim>();
        ar &binary_cell_id;
      }

    ar &data;
  }



  template <int dim, typename T>
  template <class Archive>
  void
  CellDataTransferBuffer<dim,T>::load (Archive &ar,
                                       const unsigned int /*version*/)
  {
    size_t n_cells;
    ar &n_cells;
    cell_ids.clear();
    cell_ids.reserve(n_cells);
    for (unsigned int c=0; c<n_cells; ++c)
      {
        CellId::binary_type value;
        ar &value;
        cell_ids.emplace_back(value);
      }
    ar &data;
  }



  template <typename DataType, typename MeshType>
  void
  exchange_cell_data_to_ghosts (const MeshType &mesh,
                                const std::function<boost::optional<DataType> (const typename MeshType::active_cell_iterator &)> &pack,
                                const std::function<void (const typename MeshType::active_cell_iterator &, const DataType &)> &unpack)
  {
#ifndef DEAL_II_WITH_MPI
    (void)mesh;
    (void)pack;
    (void)unpack;
    Assert(false, ExcMessage("GridTools::exchange_cell_data_to_ghosts() requires MPI."));
#else
    constexpr int dim = MeshType::dimension;
    constexpr int spacedim = MeshType::space_dimension;
    auto tria =
      static_cast<const parallel::Triangulation<dim, spacedim>*>(&mesh.get_triangulation());
    Assert (tria != nullptr,
            ExcMessage("The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));

    // map neighbor_id -> data_buffer where we accumulate the data to send
    typedef std::map<::types::subdomain_id, CellDataTransferBuffer<dim, DataType> >
    DestinationToBufferMap;
    DestinationToBufferMap destination_to_data_buffer_map;

    std::map<unsigned int, std::set<::types::subdomain_id> >
    vertices_with_ghost_neighbors = tria->compute_vertices_with_ghost_neighbors();

    for (auto cell : tria->active_cell_iterators())
      if (cell->is_locally_owned())
        {
          std::set<::types::subdomain_id> send_to;
          for (unsigned int v=0; v<GeometryInfo<dim>::vertices_per_cell; ++v)
            {
              const std::map<unsigned int, std::set<::types::subdomain_id> >::const_iterator
              neighbor_subdomains_of_vertex
                = vertices_with_ghost_neighbors.find (cell->vertex_index(v));

              if (neighbor_subdomains_of_vertex ==
                  vertices_with_ghost_neighbors.end())
                continue;

              Assert(neighbor_subdomains_of_vertex->second.size()!=0,
                     ExcInternalError());

              send_to.insert(neighbor_subdomains_of_vertex->second.begin(),
                             neighbor_subdomains_of_vertex->second.end());
            }

          if (send_to.size() > 0)
            {
              // this cell's data needs to be sent to someone
              typename MeshType::active_cell_iterator
              mesh_it (tria, cell->level(), cell->index(), &mesh);

              const boost::optional<DataType> data = pack(mesh_it);

              if (data)
                {
                  const CellId cellid = cell->id();

                  for (auto it : send_to)
                    {
                      const ::types::subdomain_id subdomain = it;

                      // find the data buffer for proc "subdomain" if it exists
                      // or create an empty one otherwise
                      typename DestinationToBufferMap::iterator p
                        = destination_to_data_buffer_map.insert (std::make_pair(subdomain,
                                                                                CellDataTransferBuffer<dim, DataType>()))
                          .first;

                      p->second.cell_ids.emplace_back(cellid);
                      p->second.data.emplace_back(data.get());
                    }
                }
            }
        }


    // 2. send our messages
    std::set<::types::subdomain_id> ghost_owners = tria->ghost_owners();
    const unsigned int n_ghost_owners = ghost_owners.size();
    std::vector<std::vector<char> > sendbuffers (n_ghost_owners);
    std::vector<MPI_Request> requests (n_ghost_owners);

    unsigned int idx=0;
    for (auto it = ghost_owners.begin();
         it!=ghost_owners.end();
         ++it, ++idx)
      {
        CellDataTransferBuffer<dim, DataType> &data = destination_to_data_buffer_map[*it];

        // pack all the data into the buffer for this recipient and send it.
        // keep data around till we can make sure that the packet has been
        // received
        sendbuffers[idx] = Utilities::pack(data);
        const int ierr = MPI_Isend(sendbuffers[idx].data(), sendbuffers[idx].size(),
                                   MPI_BYTE, *it,
                                   786, tria->get_communicator(), &requests[idx]);
        AssertThrowMPI(ierr);
      }

    // 3. receive messages
    std::vector<char> receive;
    for (unsigned int idx=0; idx<n_ghost_owners; ++idx)
      {
        MPI_Status status;
        int len;
        int ierr = MPI_Probe(MPI_ANY_SOURCE, 786, tria->get_communicator(), &status);
        AssertThrowMPI(ierr);
        ierr = MPI_Get_count(&status, MPI_BYTE, &len);
        AssertThrowMPI(ierr);

        receive.resize(len);

        char *ptr = receive.data();
        ierr = MPI_Recv(ptr, len, MPI_BYTE, status.MPI_SOURCE, status.MPI_TAG,
                        tria->get_communicator(), &status);
        AssertThrowMPI(ierr);

        auto cellinfo = Utilities::unpack<CellDataTransferBuffer<dim, DataType> >(receive);

        DataType *data = cellinfo.data.data();
        for (unsigned int c=0; c<cellinfo.cell_ids.size(); ++c, ++data)
          {
            const typename Triangulation<dim,spacedim>::cell_iterator
            tria_cell = cellinfo.cell_ids[c].to_cell(*tria);

            const typename MeshType::active_cell_iterator
            cell (tria, tria_cell->level(), tria_cell->index(), &mesh);

            unpack(cell, *data);
          }
      }

    // make sure that all communication is finished
    // when we leave this function.
    if (requests.size())
      {
        const int ierr = MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
        AssertThrowMPI(ierr);
      }
#endif // DEAL_II_WITH_MPI
  }
}

#endif

DEAL_II_NAMESPACE_CLOSE

/*----------------------------   grid_tools.h     ---------------------------*/
/* end of #ifndef dealii_grid_tools_h */
#endif
/*----------------------------   grid_tools.h     ---------------------------*/