grid_tools.h

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 2001 - 2024 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Part of the source code is dual licensed under Apache-2.0 WITH
// LLVM-exception OR LGPL-2.1-or-later. Detailed license information
// governing the source code and code contributions can be found in
// LICENSE.md and CONTRIBUTING.md at the top level directory of deal.II.
//
// ------------------------------------------------------------------------
 
#ifndef dealii_grid_tools_h
#define dealii_grid_tools_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/bounding_box.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/point.h>
 
#include <deal.II/boost_adaptors/bounding_box.h>
 
#include <deal.II/distributed/shared_tria.h>
 
#include <deal.II/dofs/dof_handler.h>
 
#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping.h>
 
#include <deal.II/grid/grid_tools_geometry.h>
#include <deal.II/grid/grid_tools_topology.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/lac/la_parallel_vector.h>
#include <deal.II/lac/petsc_vector.h>
#include <deal.II/lac/sparsity_tools.h>
#include <deal.II/lac/trilinos_vector.h>
 
#include <deal.II/numerics/rtree.h>
 
#include <boost/archive/binary_iarchive.hpp>
#include <boost/archive/binary_oarchive.hpp>
#include <boost/random/mersenne_twister.hpp>
#include <boost/serialization/array.hpp>
#include <boost/serialization/vector.hpp>
 
#ifdef DEAL_II_WITH_ZLIB
#  include <boost/iostreams/device/back_inserter.hpp>
#  include <boost/iostreams/filter/gzip.hpp>
#  include <boost/iostreams/filtering_stream.hpp>
#  include <boost/iostreams/stream.hpp>
#endif
 
#include <optional>
#include <set>
 
#ifdef DEAL_II_HAVE_CXX20
#  include <concepts>
#endif
 
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
namespace parallel
{
  namespace distributed
  {
    template <int dim, int spacedim>
    DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
    class Triangulation;
  }
} // namespace parallel
 
namespace hp
{
  template <int, int>
  class MappingCollection;
}
 
class SparsityPattern;
 
namespace GridTools
{
  template <int dim, int spacedim>
  class Cache;
}
#endif
 
namespace internal
{
  template <int dim, int spacedim, typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  class ActiveCellIterator
  {
  public:
#ifndef _MSC_VER
    using type = typename MeshType::active_cell_iterator;
#else
    using type = TriaActiveIterator<::CellAccessor<dim, spacedim>>;
#endif
  };
 
#ifdef _MSC_VER
  template <int dim, int spacedim>
  class ActiveCellIterator<dim, spacedim, DoFHandler<dim, spacedim>>
  {
  public:
    using type =
      TriaActiveIterator<::DoFCellAccessor<dim, spacedim, false>>;
  };
#endif
} // namespace internal
 
namespace GridTools
{
  template <int dim, typename Transformation, int spacedim>
  DEAL_II_CXX20_REQUIRES(
    (std::invocable<Transformation, Point<spacedim>> &&
     std::assignable_from<
       Point<spacedim> &,
       std::invoke_result_t<Transformation, Point<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);
 
 
  template <int dim, int spacedim>
  void
  rotate(const double angle, Triangulation<dim, spacedim> &triangulation);
 
  template <int dim>
  void
  rotate(const Tensor<1, 3, double> &axis,
         const double                angle,
         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>
  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,
    const unsigned int            seed = boost::random::mt19937::default_seed);
 
  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<unsigned int>>
#else
  return_type
#endif
  compute_point_locations_try_all(
    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,
    const double                                           tolerance = 1e-10,
    const std::vector<bool>                               &marked_vertices = {},
    const bool enforce_unique_mapping = true);
 
  namespace internal
  {
    template <int dim, int spacedim>
    struct DistributedComputePointLocationsInternal
    {
      DistributedComputePointLocationsInternal();
 
      void
      finalize_setup();
 
      unsigned int n_searched_points;
 
      std::vector<std::tuple<std::pair<int, int>,
                             unsigned int,
                             unsigned int,
                             Point<dim>,
                             Point<spacedim>,
                             unsigned int>>
        send_components;
 
      std::vector<unsigned int> send_ranks;
 
      std::vector<unsigned int> send_ptrs;
 
      std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
        recv_components;
 
      std::vector<unsigned int> recv_ranks;
 
      std::vector<unsigned int> recv_ptrs;
    };
 
    template <int dim, int spacedim>
    DistributedComputePointLocationsInternal<dim, spacedim>
    distributed_compute_point_locations(
      const GridTools::Cache<dim, spacedim>                 &cache,
      const std::vector<Point<spacedim>>                    &points,
      const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
      const std::vector<bool>                               &marked_vertices,
      const double                                           tolerance,
      const bool                                             perform_handshake,
      const bool enforce_unique_mapping = false);
 
 
    template <int structdim, int spacedim>
    struct DistributedComputeIntersectionLocationsInternal
    {
      using IntersectionType =
        std::array<::Point<spacedim>, structdim + 1>;
 
      std::vector<std::tuple<std::pair<int, int>,
                             unsigned int,
                             unsigned int,
                             IntersectionType>>
        send_components;
 
      std::vector<std::tuple<unsigned int, unsigned int, IntersectionType>>
        recv_components;
 
      std::vector<unsigned int> recv_ptrs;
 
      template <int dim>
      GridTools::internal::DistributedComputePointLocationsInternal<dim,
                                                                    spacedim>
      convert_to_distributed_compute_point_locations_internal(
        const unsigned int                  n_points_1D,
        const Triangulation<dim, spacedim> &tria,
        const Mapping<dim, spacedim>       &mapping,
        std::vector<Quadrature<spacedim>>  *mapped_quadratures_recv_comp =
          nullptr,
        const bool consistent_numbering_of_sender_and_receiver = false) const;
 
    private:
      std::map<unsigned int, std::vector<unsigned int>>
      communicate_indices(
        const std::vector<std::tuple<unsigned int, unsigned int, unsigned int>>
                      &point_recv_components,
        const MPI_Comm comm) const;
    };
 
    template <int structdim, int dim, int spacedim>
    DistributedComputeIntersectionLocationsInternal<structdim, spacedim>
    distributed_compute_intersection_locations(
      const Cache<dim, spacedim>                      &cache,
      const std::vector<std::vector<Point<spacedim>>> &intersection_requests,
      const std::vector<std::vector<BoundingBox<spacedim>>> &global_bboxes,
      const std::vector<bool>                               &marked_vertices,
      const double                                           tolerance);
 
  } // namespace internal
 
  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>
  DEAL_II_CXX20_REQUIRES(
    (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
  unsigned int find_closest_vertex(
    const MeshType<dim, spacedim> &mesh,
    const Point<spacedim>         &p,
    const std::vector<bool>       &marked_vertices = {});
 
  template <int dim, template <int, int> class MeshType, int spacedim>
  DEAL_II_CXX20_REQUIRES(
    (concepts::is_triangulation_or_dof_handler<MeshType<dim, 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 = {});
 
 
  template <typename MeshType>
  DEAL_II_CXX20_REQUIRES((concepts::is_triangulation_or_dof_handler<MeshType>))
#ifndef _MSC_VER
  std::vector<typename MeshType::active_cell_iterator>
#else
  std::vector<
    typename ::internal::ActiveCellIterator<MeshType::dimension,
                                                  MeshType::space_dimension,
                                                  MeshType>::type>
#endif
    find_cells_adjacent_to_vertex(const MeshType    &container,
                                  const unsigned int vertex_index);
 
  template <int dim, template <int, int> class MeshType, int spacedim>
  DEAL_II_CXX20_REQUIRES(
    (concepts::is_triangulation_or_dof_handler<MeshType<dim, 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 = {},
                                  const double             tolerance = 1.e-10);
 
  template <int dim, template <int, int> class MeshType, int spacedim>
  DEAL_II_CXX20_REQUIRES(
    (concepts::is_triangulation_or_dof_handler<MeshType<dim, 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 = {},
                                  const double             tolerance = 1.e-10);
 
  template <int dim, int spacedim>
  std::pair<typename DoFHandler<dim, spacedim>::active_cell_iterator,
            Point<dim>>
  find_active_cell_around_point(
    const hp::MappingCollection<dim, spacedim> &mapping,
    const DoFHandler<dim, spacedim>            &mesh,
    const Point<spacedim>                      &p,
    const double                                tolerance = 1.e-10);
 
  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 = {},
    const double             tolerance       = 1.e-10);
 
  template <int dim, template <int, int> class MeshType, int spacedim>
  DEAL_II_CXX20_REQUIRES(
    (concepts::is_triangulation_or_dof_handler<MeshType<dim, 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 = {},
      const RTree<std::pair<Point<spacedim>, unsigned int>> &
        used_vertices_rtree = RTree<std::pair<Point<spacedim>, unsigned int>>{},
      const double tolerance = 1.e-10,
      const RTree<
        std::pair<BoundingBox<spacedim>,
                  typename Triangulation<dim, spacedim>::active_cell_iterator>>
        *relevant_cell_bounding_boxes_rtree = nullptr);
 
  template <int dim, template <int, int> class MeshType, int spacedim>
  DEAL_II_CXX20_REQUIRES(
    (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
#ifndef _MSC_VER
  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
                        Point<dim>>>
#else
  std::vector<std::pair<
    typename ::internal::
      ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
    Point<dim>>>
#endif
    find_all_active_cells_around_point(
      const Mapping<dim, spacedim>  &mapping,
      const MeshType<dim, spacedim> &mesh,
      const Point<spacedim>         &p,
      const double                   tolerance,
      const std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
                      Point<dim>>   &first_cell,
      const std::vector<
        std::set<typename MeshType<dim, spacedim>::active_cell_iterator>>
        *vertex_to_cells = nullptr);
 
  template <int dim, template <int, int> class MeshType, int spacedim>
  DEAL_II_CXX20_REQUIRES(
    (concepts::is_triangulation_or_dof_handler<MeshType<dim, spacedim>>))
#ifndef _MSC_VER
  std::vector<std::pair<typename MeshType<dim, spacedim>::active_cell_iterator,
                        Point<dim>>>
#else
  std::vector<std::pair<
    typename ::internal::
      ActiveCellIterator<dim, spacedim, MeshType<dim, spacedim>>::type,
    Point<dim>>>
#endif
    find_all_active_cells_around_point(
      const Mapping<dim, spacedim>  &mapping,
      const MeshType<dim, spacedim> &mesh,
      const Point<spacedim>         &p,
      const double                   tolerance       = 1e-10,
      const std::vector<bool>       &marked_vertices = {});
 
  template <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  std::vector<typename MeshType::active_cell_iterator> get_active_child_cells(
    const typename MeshType::cell_iterator &cell);
 
  template <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  void get_active_neighbors(
    const typename MeshType::active_cell_iterator        &cell,
    std::vector<typename MeshType::active_cell_iterator> &active_neighbors);
 
  template <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<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 <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<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 <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  std::vector<
    typename MeshType::
      active_cell_iterator> compute_ghost_cell_halo_layer(const MeshType &mesh);
 
  template <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<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 <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  std::
    vector<typename MeshType::active_cell_iterator> compute_ghost_cell_layer_within_distance(
      const MeshType &mesh,
      const double    layer_thickness);
 
  template <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<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 spacedim>
#ifndef DOXYGEN
  std::tuple<std::map<unsigned int, 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 RTree<std::pair<BoundingBox<spacedim>, unsigned int>> &covering_rtree,
    const std::vector<Point<spacedim>>                          &points);
 
  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,
    const Mapping<dim, spacedim>                                      &mapping =
      (ReferenceCells::get_hypercube<dim>()
#ifndef _MSC_VER
         .template get_default_linear_mapping<dim, spacedim>()
#else
         .ReferenceCell::get_default_linear_mapping<dim, spacedim>()
#endif
         ));
 
  template <int dim, int spacedim>
  std::map<unsigned int, types::global_vertex_index>
  compute_local_to_global_vertex_index_map(
    const Triangulation<dim, spacedim> &triangulation);
 
  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,
                                 const bool group_siblings = true);
 
  template <int dim, int spacedim>
  void
  partition_multigrid_levels(Triangulation<dim, spacedim> &triangulation);
 
  template <int dim, int spacedim>
  std::vector<types::subdomain_id>
  get_subdomain_association(const Triangulation<dim, spacedim> &triangulation,
                            const std::vector<CellId>          &cell_ids);
 
  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 <int dim, int spacedim>
  DEAL_II_DEPRECATED typename Triangulation<dim, spacedim>::DistortedCellList
  fix_up_distorted_child_cells(
    const typename Triangulation<dim, spacedim>::DistortedCellList
                                 &distorted_cells,
    Triangulation<dim, spacedim> &triangulation);
 
 
 
  template <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<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 <int dim, int spacedim>
  std::map<
    types::global_dof_index,
    std::vector<typename DoFHandler<dim, spacedim>::active_cell_iterator>>
  get_dof_to_support_patch_map(DoFHandler<dim, spacedim> &dof_handler);
 
 
  template <typename CellIterator>
  struct PeriodicFacePair
  {
    CellIterator cell[2];
 
    unsigned int face_idx[2];
 
    unsigned char orientation;
 
    FullMatrix<double> matrix;
 
    std::size_t
    memory_consumption() const;
  };
 
 
  template <typename FaceIterator>
  std::optional<unsigned char>
  orthogonal_equality(
    const FaceIterator                                           &face1,
    const FaceIterator                                           &face2,
    const unsigned 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>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  void collect_periodic_faces(
    const MeshType          &mesh,
    const types::boundary_id b_id1,
    const types::boundary_id b_id2,
    const unsigned 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>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  void collect_periodic_faces(
    const MeshType          &mesh,
    const types::boundary_id b_id,
    const unsigned 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 <int dim, int spacedim>
  void
  assign_co_dimensional_manifold_indicators(
    Triangulation<dim, spacedim>             &tria,
    const std::function<types::manifold_id(
      const std::set<types::manifold_id> &)> &disambiguation_function =
      [](const std::set<types::manifold_id> &manifold_ids) {
        if (manifold_ids.size() == 1)
          return *manifold_ids.begin();
        else
          return numbers::flat_manifold_id;
      },
    bool overwrite_only_flat_manifold_ids = true);
  template <typename DataType, typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  void exchange_cell_data_to_ghosts(
    const MeshType                                      &mesh,
    const std::function<std::optional<DataType>(
      const typename MeshType::active_cell_iterator &)> &pack,
    const std::function<void(const typename MeshType::active_cell_iterator &,
                             const DataType &)>         &unpack,
    const std::function<bool(const typename MeshType::active_cell_iterator &)>
      &cell_filter =
        always_return<typename MeshType::active_cell_iterator, bool>{true});
 
  template <typename DataType, typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  void exchange_cell_data_to_level_ghosts(
    const MeshType                                     &mesh,
    const std::function<std::optional<DataType>(
      const typename MeshType::level_cell_iterator &)> &pack,
    const std::function<void(const typename MeshType::level_cell_iterator &,
                             const DataType &)>        &unpack,
    const std::function<bool(const typename MeshType::level_cell_iterator &)> &
      cell_filter = always_return<typename MeshType::level_cell_iterator, bool>{
        true});
 
  /* 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,
    const MPI_Comm                            mpi_communicator);
 
  template <int spacedim>
  RTree<std::pair<BoundingBox<spacedim>, unsigned int>>
  build_global_description_tree(
    const std::vector<BoundingBox<spacedim>> &local_description,
    const MPI_Comm                            mpi_communicator);
 
  template <int dim, int spacedim>
  void
  collect_coinciding_vertices(
    const Triangulation<dim, spacedim>                &tria,
    std::map<unsigned int, std::vector<unsigned int>> &coinciding_vertex_groups,
    std::map<unsigned int, unsigned int> &vertex_to_coinciding_vertex_group);
 
  template <int dim, int spacedim>
  std::map<unsigned int, std::set<::types::subdomain_id>>
  compute_vertices_with_ghost_neighbors(
    const Triangulation<dim, spacedim> &tria);
 
  template <int dim, typename VectorType>
  class MarchingCubeAlgorithm
  {
  public:
    using value_type = typename VectorType::value_type;
 
    MarchingCubeAlgorithm(const Mapping<dim, dim>       &mapping,
                          const FiniteElement<dim, dim> &fe,
                          const unsigned int             n_subdivisions = 1,
                          const double                   tolerance = 1e-10);
 
    void
    process(const DoFHandler<dim>   &background_dof_handler,
            const VectorType        &ls_vector,
            const double             iso_level,
            std::vector<Point<dim>> &vertices,
            std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
 
    void
    process(const DoFHandler<dim>   &background_dof_handler,
            const VectorType        &ls_vector,
            const double             iso_level,
            std::vector<Point<dim>> &vertices) const;
 
    void
    process_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
                 const VectorType                              &ls_vector,
                 const double                                   iso_level,
                 std::vector<Point<dim>>                       &vertices,
                 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells) const;
    void
    process_cell(const typename DoFHandler<dim>::active_cell_iterator &cell,
                 const VectorType        &ls_vector,
                 const double             iso_level,
                 std::vector<Point<dim>> &vertices) const;
 
  private:
    static Quadrature<dim>
    create_quadrature_rule(const unsigned int n_subdivisions);
 
    void
    process_cell(std::vector<value_type>                       &ls_values,
                 const std::vector<Point<dim>>                 &points,
                 const double                                   iso_level,
                 std::vector<Point<dim>>                       &vertices,
                 std::vector<CellData<dim == 1 ? 1 : dim - 1>> &cells,
                 const bool write_back_cell_data = true) const;
 
    void
    process_sub_cell(const std::vector<value_type> &,
                     const std::vector<Point<1>> &,
                     const std::vector<unsigned int> &,
                     const double,
                     std::vector<Point<1>> &,
                     std::vector<CellData<1>> &,
                     const bool) const
    {
      AssertThrow(false, ExcNotImplemented());
    }
 
    void
    process_sub_cell(const std::vector<value_type>   &ls_values,
                     const std::vector<Point<2>>     &points,
                     const std::vector<unsigned int> &mask,
                     const double                     iso_level,
                     std::vector<Point<2>>           &vertices,
                     std::vector<CellData<1>>        &cells,
                     const bool write_back_cell_data) const;
 
    void
    process_sub_cell(const std::vector<value_type>   &ls_values,
                     const std::vector<Point<3>>     &points,
                     const std::vector<unsigned int> &mask,
                     const double                     iso_level,
                     std::vector<Point<3>>           &vertices,
                     std::vector<CellData<2>>        &cells,
                     const bool write_back_cell_data) const;
 
    const unsigned int n_subdivisions;
 
    const double tolerance;
 
    mutable FEValues<dim> fe_values;
  };
 
 
 
  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
                 << '.');
 
  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>
  double
  cell_measure(
    const std::vector<Point<dim>> &all_vertices,
    const unsigned int (&indices)[GeometryInfo<dim>::vertices_per_cell])
  {
    // We forward call to the ArrayView version:
    const ArrayView<const unsigned int> view(
      indices, GeometryInfo<dim>::vertices_per_cell);
    return cell_measure(all_vertices, view);
  }
 
 
 
  // This specialization is defined here so that the general template in the
  // source file doesn't need to have further 1d overloads for the internal
  // functions it calls.
  template <>
  inline Triangulation<1, 1>::DistortedCellList
  fix_up_distorted_child_cells(const Triangulation<1, 1>::DistortedCellList &,
                               Triangulation<1, 1> &)
  {
    return {};
  }
 
 
 
  template <int dim, typename Transformation, int spacedim>
  DEAL_II_CXX20_REQUIRES(
    (std::invocable<Transformation, Point<spacedim>> &&
     std::assignable_from<
       Point<spacedim> &,
       std::invoke_result_t<Transformation, Point<spacedim>>>))
  void transform(const Transformation         &transformation,
                 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 (const unsigned int v : cell->vertex_indices())
        if (treated_vertices[cell->vertex_index(v)] == false)
          {
            // transform this vertex
            cell->vertex(v) = transformation(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 (const unsigned int face : cell->face_indices())
            if (cell->face(face)->has_children() &&
                !cell->face(face)->at_boundary())
              {
                Assert(cell->reference_cell() ==
                         ReferenceCells::get_hypercube<dim>(),
                       ExcNotImplemented());
 
                // this line has children
                cell->face(face)->child(0)->vertex(1) =
                  (cell->face(face)->vertex(0) + cell->face(face)->vertex(1)) /
                  2.0;
              }
      }
    else if (dim == 3)
      {
        typename Triangulation<dim, spacedim>::active_cell_iterator
          cell = triangulation.begin_active(),
          endc = triangulation.end();
        for (; cell != endc; ++cell)
          for (const unsigned int face : cell->face_indices())
            if (cell->face(face)->has_children() &&
                !cell->face(face)->at_boundary())
              {
                if (static_cast<uint8_t>(cell->face(face)->refinement_case()) ==
                    RefinementCase<dim - 1>::isotropic_refinement)
                  {
                    Assert(cell->reference_cell() ==
                             ReferenceCells::get_hypercube<dim>(),
                           ExcNotImplemented());
 
                    // 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;
                  }
                else
                  {
                    // Special case for anisotropic refinement
                    for (unsigned int line = 0;
                         line < GeometryInfo<dim - 1>::faces_per_cell;
                         line++)
                      if (cell->face(face)->line(line)->has_children())
                        cell->face(face)->line(line)->child(0)->vertex(1) =
                          (cell->face(face)->line(line)->vertex(0) +
                           cell->face(face)->line(line)->vertex(1)) /
                          2.0;
                  }
              }
      }
 
    // Make sure FEValues notices that the mesh has changed
    triangulation.signals.mesh_movement();
  }
 
 
 
  template <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<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 <typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<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 (const unsigned int n : cell->face_indices())
      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->is_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)->n_active_descendants();
                     ++c)
                  active_neighbors.push_back(
                    cell->neighbor_child_on_subface(n, c));
              else
                {
                  // the neighbor must be active
                  // himself
                  Assert(cell->neighbor(n)->is_active(), ExcInternalError());
                  active_neighbors.push_back(cell->neighbor(n));
                }
            }
        }
  }
 
 
 
  template <typename CellIterator>
  std::size_t
  PeriodicFacePair<CellIterator>::memory_consumption() const
  {
    return sizeof(*this) + matrix.memory_consumption();
  }
 
 
 
  namespace internal
  {
    template <typename DataType,
              typename MeshType,
              typename MeshCellIteratorType>
    DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
    inline void exchange_cell_data(
      const MeshType &mesh,
      const std::function<std::optional<DataType>(const MeshCellIteratorType &)>
        &pack,
      const std::function<void(const MeshCellIteratorType &, const DataType &)>
                                                                 &unpack,
      const std::function<bool(const MeshCellIteratorType &)>    &cell_filter,
      const std::function<void(
        const std::function<void(const MeshCellIteratorType &,
                                 const types::subdomain_id)> &)> &process_cells,
      const std::function<std::set<types::subdomain_id>(
        const parallel::TriangulationBase<MeshType::dimension,
                                          MeshType::space_dimension> &)>
        &compute_ghost_owners)
    {
#  ifndef DEAL_II_WITH_MPI
      (void)mesh;
      (void)pack;
      (void)unpack;
      (void)cell_filter;
      (void)process_cells;
      (void)compute_ghost_owners;
      Assert(false, ExcNeedsMPI());
#  else
      constexpr int dim      = MeshType::dimension;
      constexpr int spacedim = MeshType::space_dimension;
      auto          tria =
        dynamic_cast<const parallel::TriangulationBase<dim, spacedim> *>(
          &mesh.get_triangulation());
      Assert(
        tria != nullptr,
        ExcMessage(
          "The function exchange_cell_data_to_ghosts() only works with parallel triangulations."));
 
      if (const auto tria = dynamic_cast<
            const parallel::shared::Triangulation<dim, spacedim> *>(
            &mesh.get_triangulation()))
        {
          Assert(
            tria->with_artificial_cells(),
            ExcMessage(
              "The functions GridTools::exchange_cell_data_to_ghosts() and "
              "GridTools::exchange_cell_data_to_level_ghosts() can only "
              "operate on a single layer of ghost cells. However, you have "
              "given a Triangulation object of type "
              "parallel::shared::Triangulation without artificial cells "
              "resulting in an arbitrary number of ghost layers. "
              "To use this function for a Triangulation object of type "
              "parallel::shared::Triangulation, make sure to create the "
              "Triangulation object with allow_artificial_cells set to true. "
              "This results in a parallel::shared::Triangulation with only "
              "a single layer of ghost cells."));
        }
 
      // build list of cells to request for each neighbor
      std::set<::types::subdomain_id> ghost_owners =
        compute_ghost_owners(*tria);
      std::map<::types::subdomain_id,
               std::vector<typename CellId::binary_type>>
        neighbor_cell_list;
 
      for (const auto ghost_owner : ghost_owners)
        neighbor_cell_list[ghost_owner] = {};
 
      process_cells([&](const auto &cell, const auto key) -> void {
        if (cell_filter(cell))
          {
            constexpr int spacedim = MeshType::space_dimension;
            neighbor_cell_list[key].emplace_back(
              cell->id().template to_binary<spacedim>());
          }
      });
 
      Assert(ghost_owners.size() == neighbor_cell_list.size(),
             ExcInternalError());
 
 
      // Before sending & receiving, make sure we protect this section with
      // a mutex:
      static Utilities::MPI::CollectiveMutex      mutex;
      Utilities::MPI::CollectiveMutex::ScopedLock lock(
        mutex, tria->get_mpi_communicator());
 
      const int mpi_tag =
        Utilities::MPI::internal::Tags::exchange_cell_data_request;
      const int mpi_tag_reply =
        Utilities::MPI::internal::Tags::exchange_cell_data_reply;
 
      // send our requests
      std::vector<MPI_Request> requests(ghost_owners.size());
      {
        unsigned int idx = 0;
        for (const auto &it : neighbor_cell_list)
          {
            // send the data about the relevant cells
            const int ierr = MPI_Isend(it.second.data(),
                                       it.second.size() * sizeof(it.second[0]),
                                       MPI_BYTE,
                                       it.first,
                                       mpi_tag,
                                       tria->get_mpi_communicator(),
                                       &requests[idx]);
            AssertThrowMPI(ierr);
            ++idx;
          }
      }
 
      // receive requests and reply with the results
      std::vector<MPI_Request>       reply_requests(ghost_owners.size());
      std::vector<std::vector<char>> sendbuffers(ghost_owners.size());
 
      for (unsigned int idx = 0; idx < ghost_owners.size(); ++idx)
        {
          MPI_Status status;
          int        ierr = MPI_Probe(MPI_ANY_SOURCE,
                               mpi_tag,
                               tria->get_mpi_communicator(),
                               &status);
          AssertThrowMPI(ierr);
 
          int len;
          ierr = MPI_Get_count(&status, MPI_BYTE, &len);
          AssertThrowMPI(ierr);
          Assert(len % sizeof(typename CellId::binary_type) == 0,
                 ExcInternalError());
 
          const unsigned int n_cells =
            len / sizeof(typename CellId::binary_type);
          std::vector<typename CellId::binary_type> cells_with_requests(
            n_cells);
          std::vector<DataType> data_to_send;
          data_to_send.reserve(n_cells);
          std::vector<bool> cell_carries_data(n_cells, false);
 
          ierr = MPI_Recv(cells_with_requests.data(),
                          len,
                          MPI_BYTE,
                          status.MPI_SOURCE,
                          status.MPI_TAG,
                          tria->get_mpi_communicator(),
                          &status);
          AssertThrowMPI(ierr);
 
          // store data for each cell
          for (unsigned int c = 0; c < static_cast<unsigned int>(n_cells); ++c)
            {
              const auto cell =
                tria->create_cell_iterator(CellId(cells_with_requests[c]));
 
              MeshCellIteratorType mesh_it(tria,
                                           cell->level(),
                                           cell->index(),
                                           &mesh);
 
              const std::optional<DataType> data = pack(mesh_it);
              if (data)
                {
                  data_to_send.emplace_back(*data);
                  cell_carries_data[c] = true;
                }
            }
 
          // collect data for sending the reply in a buffer
 
          // (a) make room for storing the local offsets in case we receive
          // other data
          sendbuffers[idx].resize(sizeof(std::size_t));
 
          // (b) append the actual data and store how much memory it
          // corresponds to, which we then insert into the leading position of
          // the sendbuffer
          std::size_t size_of_send =
            Utilities::pack(data_to_send,
                            sendbuffers[idx],
                            /*enable_compression*/ false);
          std::memcpy(sendbuffers[idx].data(),
                      &size_of_send,
                      sizeof(std::size_t));
 
          // (c) append information of certain cells that got left out in case
          // we need it
          if (data_to_send.size() < n_cells)
            Utilities::pack(cell_carries_data,
                            sendbuffers[idx],
                            /*enable_compression*/ false);
 
          // send data
          ierr = MPI_Isend(sendbuffers[idx].data(),
                           sendbuffers[idx].size(),
                           MPI_BYTE,
                           status.MPI_SOURCE,
                           mpi_tag_reply,
                           tria->get_mpi_communicator(),
                           &reply_requests[idx]);
          AssertThrowMPI(ierr);
        }
 
      // finally receive the replies
      std::vector<char> receive;
      for (unsigned int id = 0; id < neighbor_cell_list.size(); ++id)
        {
          MPI_Status status;
          int        ierr = MPI_Probe(MPI_ANY_SOURCE,
                               mpi_tag_reply,
                               tria->get_mpi_communicator(),
                               &status);
          AssertThrowMPI(ierr);
 
          int len;
          ierr = MPI_Get_count(&status, MPI_BYTE, &len);
          AssertThrowMPI(ierr);
 
          receive.resize(len);
 
          ierr = MPI_Recv(receive.data(),
                          len,
                          MPI_BYTE,
                          status.MPI_SOURCE,
                          status.MPI_TAG,
                          tria->get_mpi_communicator(),
                          &status);
          AssertThrowMPI(ierr);
 
          // (a) first determine the length of the data section in the
          // received buffer
          auto        data_iterator = receive.begin();
          std::size_t size_of_received_data =
            Utilities::unpack<std::size_t>(data_iterator,
                                           data_iterator + sizeof(std::size_t));
          data_iterator += sizeof(std::size_t);
 
          // (b) unpack the data section in the indicated region
          auto received_data = Utilities::unpack<std::vector<DataType>>(
            data_iterator,
            data_iterator + size_of_received_data,
            /*enable_compression*/ false);
          data_iterator += size_of_received_data;
 
          // (c) check if the received data contained fewer entries than the
          // number of cells we identified in the beginning, in which case we
          // need to extract the boolean vector with the relevant information
          const std::vector<typename CellId::binary_type> &this_cell_list =
            neighbor_cell_list[status.MPI_SOURCE];
          AssertIndexRange(received_data.size(), this_cell_list.size() + 1);
          std::vector<bool> cells_with_data;
          if (received_data.size() < this_cell_list.size())
            {
              cells_with_data = Utilities::unpack<std::vector<bool>>(
                data_iterator, receive.end(), /*enable_compression*/ false);
              AssertDimension(cells_with_data.size(), this_cell_list.size());
            }
 
          // (d) go through the received data and call the user-provided
          // unpack function
          auto received_data_iterator = received_data.begin();
          for (unsigned int c = 0; c < this_cell_list.size(); ++c)
            if (cells_with_data.empty() || cells_with_data[c])
              {
                const typename Triangulation<dim, spacedim>::cell_iterator
                  tria_cell = tria->create_cell_iterator(this_cell_list[c]);
 
                MeshCellIteratorType cell(tria,
                                          tria_cell->level(),
                                          tria_cell->index(),
                                          &mesh);
 
                unpack(cell, *received_data_iterator);
                ++received_data_iterator;
              }
        }
 
      // make sure that all communication is finished
      // when we leave this function.
      if (requests.size() > 0)
        {
          const int ierr =
            MPI_Waitall(requests.size(), requests.data(), MPI_STATUSES_IGNORE);
          AssertThrowMPI(ierr);
        }
      if (reply_requests.size() > 0)
        {
          const int ierr = MPI_Waitall(reply_requests.size(),
                                       reply_requests.data(),
                                       MPI_STATUSES_IGNORE);
          AssertThrowMPI(ierr);
        }
 
 
#  endif // DEAL_II_WITH_MPI
    }
 
  } // namespace internal
 
  template <typename DataType, typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  inline void exchange_cell_data_to_ghosts(
    const MeshType                                      &mesh,
    const std::function<std::optional<DataType>(
      const typename MeshType::active_cell_iterator &)> &pack,
    const std::function<void(const typename MeshType::active_cell_iterator &,
                             const DataType &)>         &unpack,
    const std::function<bool(const typename MeshType::active_cell_iterator &)>
      &cell_filter)
  {
#  ifndef DEAL_II_WITH_MPI
    (void)mesh;
    (void)pack;
    (void)unpack;
    (void)cell_filter;
    Assert(false, ExcNeedsMPI());
#  else
    internal::exchange_cell_data<DataType,
                                 MeshType,
                                 typename MeshType::active_cell_iterator>(
      mesh,
      pack,
      unpack,
      cell_filter,
      [&](const auto &process) {
        for (const auto &cell : mesh.active_cell_iterators())
          if (cell->is_ghost())
            process(cell, cell->subdomain_id());
      },
      [](const auto &tria) { return tria.ghost_owners(); });
#  endif
  }
 
 
 
  template <typename DataType, typename MeshType>
  DEAL_II_CXX20_REQUIRES(concepts::is_triangulation_or_dof_handler<MeshType>)
  inline void exchange_cell_data_to_level_ghosts(
    const MeshType                                     &mesh,
    const std::function<std::optional<DataType>(
      const typename MeshType::level_cell_iterator &)> &pack,
    const std::function<void(const typename MeshType::level_cell_iterator &,
                             const DataType &)>        &unpack,
    const std::function<bool(const typename MeshType::level_cell_iterator &)>
      &cell_filter)
  {
#  ifndef DEAL_II_WITH_MPI
    (void)mesh;
    (void)pack;
    (void)unpack;
    (void)cell_filter;
    Assert(false, ExcNeedsMPI());
#  else
    internal::exchange_cell_data<DataType,
                                 MeshType,
                                 typename MeshType::level_cell_iterator>(
      mesh,
      pack,
      unpack,
      cell_filter,
      [&](const auto &process) {
        for (const auto &cell : mesh.cell_iterators())
          if (cell->is_ghost_on_level())
            process(cell, cell->level_subdomain_id());
      },
      [](const auto &tria) { return tria.level_ghost_owners(); });
#  endif
  }
} // namespace GridTools
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif