tria.cc

Go to the documentation of this file.

// ------------------------------------------------------------------------
//
// SPDX-License-Identifier: LGPL-2.1-or-later
// Copyright (C) 1999 - 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.
//
// ------------------------------------------------------------------------
 
 
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/memory_consumption.h>
#include <deal.II/base/mpi.templates.h>
#include <deal.II/base/mpi_large_count.h>
#include <deal.II/base/mpi_stub.h>
#include <deal.II/base/thread_management.h>
#include <deal.II/base/utilities.h>
 
#include <deal.II/distributed/tria.h>
 
#include <deal.II/grid/connectivity.h>
#include <deal.II/grid/grid_tools.h>
#include <deal.II/grid/magic_numbers.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_faces.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_levels.h>
 
#include <boost/archive/text_iarchive.hpp>
#include <boost/archive/text_oarchive.hpp>
 
#include <algorithm>
#include <array>
#include <cmath>
#include <cstdint>
#include <fstream>
#include <functional>
#include <limits>
#include <list>
#include <map>
#include <memory>
#include <numeric>
 
 
DEAL_II_NAMESPACE_OPEN
 
 
namespace internal
{
  namespace TriangulationImplementation
  {
    NumberCache<1>::NumberCache()
      : n_levels(0)
      , n_lines(0)
      , n_active_lines(0)
    // all other fields are
    // default constructed
    {}
 
 
 
    std::size_t
    NumberCache<1>::memory_consumption() const
    {
      std::size_t mem =
        MemoryConsumption::memory_consumption(n_levels) +
        MemoryConsumption::memory_consumption(n_lines) +
        MemoryConsumption::memory_consumption(n_lines_level) +
        MemoryConsumption::memory_consumption(n_active_lines) +
        MemoryConsumption::memory_consumption(n_active_lines_level);
 
      if (active_cell_index_partitioner)
        mem += active_cell_index_partitioner->memory_consumption();
 
      for (const auto &partitioner : level_cell_index_partitioners)
        if (partitioner)
          mem += partitioner->memory_consumption();
 
      return mem;
    }
 
 
    NumberCache<2>::NumberCache()
      : n_quads(0)
      , n_active_quads(0)
    // all other fields are
    // default constructed
    {}
 
 
 
    std::size_t
    NumberCache<2>::memory_consumption() const
    {
      return (NumberCache<1>::memory_consumption() +
              MemoryConsumption::memory_consumption(n_quads) +
              MemoryConsumption::memory_consumption(n_quads_level) +
              MemoryConsumption::memory_consumption(n_active_quads) +
              MemoryConsumption::memory_consumption(n_active_quads_level));
    }
 
 
 
    NumberCache<3>::NumberCache()
      : n_hexes(0)
      , n_active_hexes(0)
    // all other fields are
    // default constructed
    {}
 
 
 
    std::size_t
    NumberCache<3>::memory_consumption() const
    {
      return (NumberCache<2>::memory_consumption() +
              MemoryConsumption::memory_consumption(n_hexes) +
              MemoryConsumption::memory_consumption(n_hexes_level) +
              MemoryConsumption::memory_consumption(n_active_hexes) +
              MemoryConsumption::memory_consumption(n_active_hexes_level));
    }
  } // namespace TriangulationImplementation
 
 
  template <int dim, int spacedim>
  DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
  CellAttachedDataSerializer<dim, spacedim>::CellAttachedDataSerializer()
    : variable_size_data_stored(false)
  {}
 
 
  template <int dim, int spacedim>
  DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
  void CellAttachedDataSerializer<dim, spacedim>::pack_data(
    const std::vector<cell_relation_t> &cell_relations,
    const std::vector<
      typename internal::CellAttachedData<dim, spacedim>::pack_callback_t>
      &pack_callbacks_fixed,
    const std::vector<
      typename internal::CellAttachedData<dim, spacedim>::pack_callback_t>
                   &pack_callbacks_variable,
    const MPI_Comm &mpi_communicator)
  {
    Assert(src_data_fixed.empty(),
           ExcMessage("Previously packed data has not been released yet!"));
    Assert(src_sizes_variable.empty(), ExcInternalError());
 
    const unsigned int n_callbacks_fixed    = pack_callbacks_fixed.size();
    const unsigned int n_callbacks_variable = pack_callbacks_variable.size();
 
    // Store information that we packed variable size data in
    // a member variable for later.
    variable_size_data_stored = (n_callbacks_variable > 0);
 
    // If variable transfer is scheduled, we will store the data size that
    // each variable size callback function writes in this auxiliary
    // container. The information will be stored by each cell in this vector
    // temporarily.
    std::vector<unsigned int> cell_sizes_variable_cumulative(
      n_callbacks_variable);
 
    // Prepare the buffer structure, in which each callback function will
    // store its data for each active cell.
    // The outmost shell in this container construct corresponds to the
    // data packed per cell. The next layer resembles the data that
    // each callback function packs on the corresponding cell. These
    // buffers are chains of chars stored in an std::vector<char>.
    // A visualisation of the data structure:
    /* clang-format off */
      // |             cell_1                | |             cell_2                | ...
      // ||  callback_1  ||  callback_2  |...| ||  callback_1  ||  callback_2  |...| ...
      // |||char|char|...|||char|char|...|...| |||char|char|...|||char|char|...|...| ...
    /* clang-format on */
    std::vector<std::vector<std::vector<char>>> packed_fixed_size_data(
      cell_relations.size());
    std::vector<std::vector<std::vector<char>>> packed_variable_size_data(
      variable_size_data_stored ? cell_relations.size() : 0);
 
    //
    // --------- Pack data for fixed and variable size transfer ---------
    //
    // Iterate over all cells, call all callback functions on each cell,
    // and store their data in the corresponding buffer scope.
    {
      auto cell_rel_it           = cell_relations.cbegin();
      auto data_cell_fixed_it    = packed_fixed_size_data.begin();
      auto data_cell_variable_it = packed_variable_size_data.begin();
      for (; cell_rel_it != cell_relations.cend(); ++cell_rel_it)
        {
          const auto &dealii_cell = cell_rel_it->first;
          const auto &cell_status = cell_rel_it->second;
 
          // Assertions about the tree structure.
          switch (cell_status)
            {
              case CellStatus::cell_will_persist:
              case CellStatus::cell_will_be_refined:
                // double check the condition that we will only ever attach
                // data to active cells when we get here
                Assert(dealii_cell->is_active(), ExcInternalError());
                break;
 
              case CellStatus::children_will_be_coarsened:
                // double check the condition that we will only ever attach
                // data to cells with children when we get here. however, we
                // can only tolerate one level of coarsening at a time, so
                // check that the children are all active
                Assert(dealii_cell->is_active() == false, ExcInternalError());
                for (unsigned int c = 0;
                     c < GeometryInfo<dim>::max_children_per_cell;
                     ++c)
                  Assert(dealii_cell->child(c)->is_active(),
                         ExcInternalError());
                break;
 
              case CellStatus::cell_invalid:
                // do nothing on invalid cells
                break;
 
              default:
                DEAL_II_ASSERT_UNREACHABLE();
                break;
            }
 
          // Reserve memory corresponding to the number of callback
          // functions that will be called.
          // If variable size transfer is scheduled, we need to leave
          // room for an array that holds information about how many
          // bytes each of the variable size callback functions will
          // write.
          // On cells flagged with CellStatus::cell_invalid, only its CellStatus
          // will be stored.
          const unsigned int n_fixed_size_data_sets_on_cell =
            1 + ((cell_status == CellStatus::cell_invalid) ?
                   0 :
                   ((variable_size_data_stored ? 1 : 0) + n_callbacks_fixed));
          data_cell_fixed_it->resize(n_fixed_size_data_sets_on_cell);
 
          // We continue with packing all data on this specific cell.
          auto data_fixed_it = data_cell_fixed_it->begin();
 
          // First, we pack the CellStatus information.
          // to get consistent data sizes on each cell for the fixed size
          // transfer, we won't allow compression
          *data_fixed_it =
            Utilities::pack(cell_status, /*allow_compression=*/false);
          ++data_fixed_it;
 
          // Proceed with all registered callback functions.
          // Skip cells with the CellStatus::cell_invalid flag.
          if (cell_status != CellStatus::cell_invalid)
            {
              // Pack fixed size data.
              for (auto callback_it = pack_callbacks_fixed.cbegin();
                   callback_it != pack_callbacks_fixed.cend();
                   ++callback_it, ++data_fixed_it)
                {
                  *data_fixed_it = (*callback_it)(dealii_cell, cell_status);
                }
 
              // Pack variable size data.
              // If we store variable size data, we need to transfer
              // the sizes of each corresponding callback function
              // via fixed size transfer as well.
              if (variable_size_data_stored)
                {
                  const unsigned int n_variable_size_data_sets_on_cell =
                    ((cell_status == CellStatus::cell_invalid) ?
                       0 :
                       n_callbacks_variable);
                  data_cell_variable_it->resize(
                    n_variable_size_data_sets_on_cell);
 
                  auto callback_it      = pack_callbacks_variable.cbegin();
                  auto data_variable_it = data_cell_variable_it->begin();
                  auto sizes_variable_it =
                    cell_sizes_variable_cumulative.begin();
                  for (; callback_it != pack_callbacks_variable.cend();
                       ++callback_it, ++data_variable_it, ++sizes_variable_it)
                    {
                      *data_variable_it =
                        (*callback_it)(dealii_cell, cell_status);
 
                      // Store data sizes for each callback function first.
                      // Make it cumulative below.
                      *sizes_variable_it = data_variable_it->size();
                    }
 
                  // Turn size vector into its cumulative representation.
                  std::partial_sum(cell_sizes_variable_cumulative.begin(),
                                   cell_sizes_variable_cumulative.end(),
                                   cell_sizes_variable_cumulative.begin());
 
                  // Serialize cumulative variable size vector
                  // value-by-value. This way we can circumvent the overhead
                  // of storing the container object as a whole, since we
                  // know its size by the number of registered callback
                  // functions.
                  data_fixed_it->resize(n_callbacks_variable *
                                        sizeof(unsigned int));
                  for (unsigned int i = 0; i < n_callbacks_variable; ++i)
                    std::memcpy(&(data_fixed_it->at(i * sizeof(unsigned int))),
                                &(cell_sizes_variable_cumulative.at(i)),
                                sizeof(unsigned int));
 
                  ++data_fixed_it;
                }
 
              // Double check that we packed everything we wanted
              // in the fixed size buffers.
              Assert(data_fixed_it == data_cell_fixed_it->end(),
                     ExcInternalError());
            }
 
          ++data_cell_fixed_it;
 
          // Increment the variable size data iterator
          // only if we actually pack this kind of data
          // to avoid getting out of bounds.
          if (variable_size_data_stored)
            ++data_cell_variable_it;
        } // loop over cell_relations
    }
 
    //
    // ----------- Gather data sizes for fixed size transfer ------------
    //
    // Generate a vector which stores the sizes of each callback function,
    // including the packed CellStatus transfer.
    // Find the very first cell that we wrote to with all callback
    // functions (i.e. a cell that was not flagged with
    // CellStatus::cell_invalid) and store the sizes of each buffer.
    //
    // To deal with the case that at least one of the processors does not
    // own any cell at all, we will exchange the information about the data
    // sizes among them later. The code in between is still well-defined,
    // since the following loops will be skipped.
    std::vector<unsigned int> local_sizes_fixed(
      1 + n_callbacks_fixed + (variable_size_data_stored ? 1 : 0));
    for (const auto &data_cell : packed_fixed_size_data)
      {
        if (data_cell.size() == local_sizes_fixed.size())
          {
            auto sizes_fixed_it = local_sizes_fixed.begin();
            auto data_fixed_it  = data_cell.cbegin();
            for (; data_fixed_it != data_cell.cend();
                 ++data_fixed_it, ++sizes_fixed_it)
              {
                *sizes_fixed_it = data_fixed_it->size();
              }
 
            break;
          }
      }
 
    // Check if all cells have valid sizes.
    for (auto data_cell_fixed_it = packed_fixed_size_data.cbegin();
         data_cell_fixed_it != packed_fixed_size_data.cend();
         ++data_cell_fixed_it)
      {
        Assert((data_cell_fixed_it->size() == 1) ||
                 (data_cell_fixed_it->size() == local_sizes_fixed.size()),
               ExcInternalError());
      }
 
    // Share information about the packed data sizes
    // of all callback functions across all processors, in case one
    // of them does not own any cells at all.
    std::vector<unsigned int> global_sizes_fixed(local_sizes_fixed.size());
    Utilities::MPI::max(local_sizes_fixed,
                        mpi_communicator,
                        global_sizes_fixed);
 
    // Construct cumulative sizes, since this is the only information
    // we need from now on.
    sizes_fixed_cumulative.resize(global_sizes_fixed.size());
    std::partial_sum(global_sizes_fixed.begin(),
                     global_sizes_fixed.end(),
                     sizes_fixed_cumulative.begin());
 
    //
    // ---------- Gather data sizes for variable size transfer ----------
    //
    if (variable_size_data_stored)
      {
        src_sizes_variable.reserve(packed_variable_size_data.size());
        for (const auto &data_cell : packed_variable_size_data)
          {
            int variable_data_size_on_cell = 0;
 
            for (const auto &data : data_cell)
              variable_data_size_on_cell += data.size();
 
            src_sizes_variable.push_back(variable_data_size_on_cell);
          }
      }
 
    //
    // ------------------------ Build buffers ---------------------------
    //
    const unsigned int expected_size_fixed =
      cell_relations.size() * sizes_fixed_cumulative.back();
    const unsigned int expected_size_variable =
      std::accumulate(src_sizes_variable.begin(),
                      src_sizes_variable.end(),
                      std::vector<int>::size_type(0));
 
    // Move every piece of packed fixed size data into the consecutive
    // buffer.
    src_data_fixed.reserve(expected_size_fixed);
    for (const auto &data_cell_fixed : packed_fixed_size_data)
      {
        // Move every fraction of packed data into the buffer
        // reserved for this particular cell.
        for (const auto &data_fixed : data_cell_fixed)
          std::move(data_fixed.begin(),
                    data_fixed.end(),
                    std::back_inserter(src_data_fixed));
 
        // If we only packed the CellStatus information
        // (i.e. encountered a cell flagged CellStatus::cell_invalid),
        // fill the remaining space with invalid entries.
        // We can skip this if there is nothing else to pack.
        if ((data_cell_fixed.size() == 1) &&
            (sizes_fixed_cumulative.size() > 1))
          {
            const std::size_t bytes_skipped =
              sizes_fixed_cumulative.back() - sizes_fixed_cumulative.front();
 
            src_data_fixed.insert(src_data_fixed.end(),
                                  bytes_skipped,
                                  static_cast<char>(-1)); // invalid_char
          }
      }
 
    // Move every piece of packed variable size data into the consecutive
    // buffer.
    if (variable_size_data_stored)
      {
        src_data_variable.reserve(expected_size_variable);
        for (const auto &data_cell : packed_variable_size_data)
          {
            // Move every fraction of packed data into the buffer
            // reserved for this particular cell.
            for (const auto &data : data_cell)
              std::move(data.begin(),
                        data.end(),
                        std::back_inserter(src_data_variable));
          }
      }
 
    // Double check that we packed everything correctly.
    Assert(src_data_fixed.size() == expected_size_fixed, ExcInternalError());
    Assert(src_data_variable.size() == expected_size_variable,
           ExcInternalError());
  }
 
 
 
  template <int dim, int spacedim>
  DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
  void CellAttachedDataSerializer<dim, spacedim>::unpack_cell_status(
    std::vector<
      typename CellAttachedDataSerializer<dim, spacedim>::cell_relation_t>
      &cell_relations) const
  {
    Assert(sizes_fixed_cumulative.size() > 0,
           ExcMessage("No data has been packed!"));
    if (cell_relations.size() > 0)
      {
        Assert(dest_data_fixed.size() > 0,
               ExcMessage("No data has been received!"));
      }
 
    // Size of CellStatus object that will be unpacked on each cell.
    const unsigned int size = sizes_fixed_cumulative.front();
 
    // Iterate over all cells and overwrite the CellStatus
    // information from the transferred data.
    // Proceed buffer iterator position to next cell after
    // each iteration.
    auto cell_rel_it   = cell_relations.begin();
    auto dest_fixed_it = dest_data_fixed.cbegin();
    for (; cell_rel_it != cell_relations.end();
         ++cell_rel_it, dest_fixed_it += sizes_fixed_cumulative.back())
      {
        cell_rel_it->second = // cell_status
          Utilities::unpack<CellStatus>(dest_fixed_it,
                                        dest_fixed_it + size,
                                        /*allow_compression=*/false);
      }
  }
 
 
 
  template <int dim, int spacedim>
  DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
  void CellAttachedDataSerializer<dim, spacedim>::unpack_data(
    const std::vector<
      typename CellAttachedDataSerializer<dim, spacedim>::cell_relation_t>
                      &cell_relations,
    const unsigned int handle,
    const std::function<
      void(const cell_iterator &,
           const CellStatus &,
           const boost::iterator_range<std::vector<char>::const_iterator> &)>
      &unpack_callback) const
  {
    // We decode the handle returned by register_data_attach() back into
    // a format we can use. All even handles belong to those callback
    // functions which write/read variable size data, all odd handles
    // interact with fixed size buffers.
    const bool         callback_variable_transfer = (handle % 2 == 0);
    const unsigned int callback_index             = handle / 2;
 
    // Cells will always receive fixed size data (i.e., CellStatus
    // information), but not necessarily variable size data (e.g., with a
    // ParticleHandler a cell might not contain any particle at all).
    // Thus it is sufficient to check if fixed size data has been received.
    Assert(sizes_fixed_cumulative.size() > 0,
           ExcMessage("No data has been packed!"));
    if (cell_relations.size() > 0)
      {
        Assert(dest_data_fixed.size() > 0,
               ExcMessage("No data has been received!"));
      }
 
    std::vector<char>::const_iterator dest_data_it;
    std::vector<char>::const_iterator dest_sizes_cell_it;
 
    // Depending on whether our callback function unpacks fixed or
    // variable size data, we have to pursue different approaches
    // to localize the correct fraction of the buffer from which
    // we are allowed to read.
    unsigned int offset         = numbers::invalid_unsigned_int;
    unsigned int size           = numbers::invalid_unsigned_int;
    unsigned int data_increment = numbers::invalid_unsigned_int;
 
    if (callback_variable_transfer)
      {
        // For the variable size data, we need to extract the
        // data size from the fixed size buffer on each cell.
        //
        // We packed this information last, so the last packed
        // object in the fixed size buffer corresponds to the
        // variable data sizes.
        //
        // The last entry of sizes_fixed_cumulative corresponds
        // to the size of all fixed size data packed on the cell.
        // To get the offset for the last packed object, we need
        // to get the next-to-last entry.
        const unsigned int offset_variable_data_sizes =
          sizes_fixed_cumulative[sizes_fixed_cumulative.size() - 2];
 
        // This iterator points to the data size that the
        // callback_function packed for each specific cell.
        // Adjust buffer iterator to the offset of the callback
        // function so that we only have to advance its position
        // to the next cell after each iteration.
        dest_sizes_cell_it = dest_data_fixed.cbegin() +
                             offset_variable_data_sizes +
                             callback_index * sizeof(unsigned int);
 
        // Let the data iterator point to the correct buffer.
        dest_data_it = dest_data_variable.cbegin();
      }
    else
      {
        // For the fixed size data, we can get the information about
        // the buffer location on each cell directly from the
        // sizes_fixed_cumulative vector.
        offset         = sizes_fixed_cumulative[callback_index];
        size           = sizes_fixed_cumulative[callback_index + 1] - offset;
        data_increment = sizes_fixed_cumulative.back();
 
        // Let the data iterator point to the correct buffer.
        // Adjust buffer iterator to the offset of the callback
        // function so that we only have to advance its position
        // to the next cell after each iteration.
        if (cell_relations.begin() != cell_relations.end())
          dest_data_it = dest_data_fixed.cbegin() + offset;
      }
 
    // Iterate over all cells and unpack the transferred data.
    auto cell_rel_it   = cell_relations.begin();
    auto dest_sizes_it = dest_sizes_variable.cbegin();
    for (; cell_rel_it != cell_relations.end(); ++cell_rel_it)
      {
        const auto &dealii_cell = cell_rel_it->first;
        const auto &cell_status = cell_rel_it->second;
 
        if (callback_variable_transfer)
          {
            // Update the increment according to the whole data size
            // of the current cell.
            data_increment = *dest_sizes_it;
 
            if (cell_status != CellStatus::cell_invalid)
              {
                // Extract the corresponding values for offset and size from
                // the cumulative sizes array stored in the fixed size
                // buffer.
                if (callback_index == 0)
                  offset = 0;
                else
                  std::memcpy(&offset,
                              &(*(dest_sizes_cell_it - sizeof(unsigned int))),
                              sizeof(unsigned int));
 
                std::memcpy(&size,
                            &(*dest_sizes_cell_it),
                            sizeof(unsigned int));
 
                size -= offset;
 
                // Move the data iterator to the corresponding position
                // of the callback function and adjust the increment
                // accordingly.
                dest_data_it += offset;
                data_increment -= offset;
              }
 
            // Advance data size iterators to the next cell, avoid iterating
            // past the end of dest_sizes_cell_it
            if (cell_rel_it != cell_relations.end() - 1)
              dest_sizes_cell_it += sizes_fixed_cumulative.back();
            ++dest_sizes_it;
          }
 
        switch (cell_status)
          {
            case CellStatus::cell_will_persist:
            case CellStatus::children_will_be_coarsened:
              unpack_callback(dealii_cell,
                              cell_status,
                              boost::make_iterator_range(dest_data_it,
                                                         dest_data_it + size));
              break;
 
            case CellStatus::cell_will_be_refined:
              unpack_callback(dealii_cell->parent(),
                              cell_status,
                              boost::make_iterator_range(dest_data_it,
                                                         dest_data_it + size));
              break;
 
            case CellStatus::cell_invalid:
              // Skip this cell.
              break;
 
            default:
              DEAL_II_ASSERT_UNREACHABLE();
              break;
          }
 
        if (cell_rel_it != cell_relations.end() - 1)
          dest_data_it += data_increment;
      }
  }
 
 
 
  template <int dim, int spacedim>
  DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
  void CellAttachedDataSerializer<dim, spacedim>::save(
    const unsigned int global_first_cell,
    const unsigned int global_num_cells,
    const std::string &file_basename,
    const MPI_Comm    &mpi_communicator) const
  {
    Assert(sizes_fixed_cumulative.size() > 0,
           ExcMessage("No data has been packed!"));
 
#ifdef DEAL_II_WITH_MPI
    // Large fractions of this function have been copied from
    // DataOutInterface::write_vtu_in_parallel.
    // TODO: Write general MPIIO interface.
 
    const unsigned int myrank =
      Utilities::MPI::this_mpi_process(mpi_communicator);
    const unsigned int mpisize =
      Utilities::MPI::n_mpi_processes(mpi_communicator);
 
    if (mpisize > 1)
      {
        const unsigned int bytes_per_cell = sizes_fixed_cumulative.back();
 
        //
        // ---------- Fixed size data ----------
        //
        {
          const std::string fname_fixed =
            std::string(file_basename) + "_fixed.data";
 
          MPI_Info info;
          int      ierr = MPI_Info_create(&info);
          AssertThrowMPI(ierr);
 
          MPI_File fh;
          ierr = MPI_File_open(mpi_communicator,
                               fname_fixed.c_str(),
                               MPI_MODE_CREATE | MPI_MODE_WRONLY,
                               info,
                               &fh);
          AssertThrowMPI(ierr);
 
          ierr = MPI_File_set_size(fh, 0); // delete the file contents
          AssertThrowMPI(ierr);
          // this barrier is necessary, because otherwise others might already
          // write while one core is still setting the size to zero.
          ierr = MPI_Barrier(mpi_communicator);
          AssertThrowMPI(ierr);
          ierr = MPI_Info_free(&info);
          AssertThrowMPI(ierr);
          // ------------------
 
          // Write cumulative sizes to file.
          // Since each processor owns the same information about the data
          // sizes, it is sufficient to let only the first processor perform
          // this task.
          if (myrank == 0)
            {
              ierr = Utilities::MPI::LargeCount::File_write_at_c(
                fh,
                0,
                sizes_fixed_cumulative.data(),
                sizes_fixed_cumulative.size(),
                MPI_UNSIGNED,
                MPI_STATUS_IGNORE);
              AssertThrowMPI(ierr);
            }
 
          // Write packed data to file simultaneously.
          const MPI_Offset size_header =
            sizes_fixed_cumulative.size() * sizeof(unsigned int);
 
          // Make sure we do the following computation in 64bit integers to be
          // able to handle 4GB+ files:
          const MPI_Offset my_global_file_position =
            size_header +
            static_cast<MPI_Offset>(global_first_cell) * bytes_per_cell;
 
          ierr =
            Utilities::MPI::LargeCount::File_write_at_c(fh,
                                                        my_global_file_position,
                                                        src_data_fixed.data(),
                                                        src_data_fixed.size(),
                                                        MPI_BYTE,
                                                        MPI_STATUS_IGNORE);
          AssertThrowMPI(ierr);
 
          ierr = MPI_File_close(&fh);
          AssertThrowMPI(ierr);
        }
 
 
 
        //
        // ---------- Variable size data ----------
        //
        if (variable_size_data_stored)
          {
            const std::string fname_variable =
              std::string(file_basename) + "_variable.data";
 
            MPI_Info info;
            int      ierr = MPI_Info_create(&info);
            AssertThrowMPI(ierr);
 
            MPI_File fh;
            ierr = MPI_File_open(mpi_communicator,
                                 fname_variable.c_str(),
                                 MPI_MODE_CREATE | MPI_MODE_WRONLY,
                                 info,
                                 &fh);
            AssertThrowMPI(ierr);
 
            ierr = MPI_File_set_size(fh, 0); // delete the file contents
            AssertThrowMPI(ierr);
            // this barrier is necessary, because otherwise others might already
            // write while one core is still setting the size to zero.
            ierr = MPI_Barrier(mpi_communicator);
            AssertThrowMPI(ierr);
            ierr = MPI_Info_free(&info);
            AssertThrowMPI(ierr);
 
            // Write sizes of each cell into file simultaneously.
            {
              const MPI_Offset my_global_file_position =
                static_cast<MPI_Offset>(global_first_cell) *
                sizeof(unsigned int);
 
              // It is very unlikely that a single process has more than
              // 2 billion cells, but we might as well check.
              AssertThrow(src_sizes_variable.size() <
                            static_cast<std::size_t>(
                              std::numeric_limits<int>::max()),
                          ExcNotImplemented());
 
              ierr = Utilities::MPI::LargeCount::File_write_at_c(
                fh,
                my_global_file_position,
                src_sizes_variable.data(),
                src_sizes_variable.size(),
                MPI_INT,
                MPI_STATUS_IGNORE);
              AssertThrowMPI(ierr);
            }
 
            // Gather size of data in bytes we want to store from this
            // processor and compute the prefix sum. We do this in 64 bit
            // to avoid overflow for files larger than 4GB:
            const std::uint64_t size_on_proc = src_data_variable.size();
            std::uint64_t       prefix_sum   = 0;
            ierr                             = MPI_Exscan(&size_on_proc,
                              &prefix_sum,
                              1,
                              MPI_UINT64_T,
                              MPI_SUM,
                              mpi_communicator);
            AssertThrowMPI(ierr);
 
            const MPI_Offset my_global_file_position =
              static_cast<MPI_Offset>(global_num_cells) * sizeof(unsigned int) +
              prefix_sum;
 
            // Write data consecutively into file.
            ierr = Utilities::MPI::LargeCount::File_write_at_c(
              fh,
              my_global_file_position,
              src_data_variable.data(),
              src_data_variable.size(),
              MPI_BYTE,
              MPI_STATUS_IGNORE);
            AssertThrowMPI(ierr);
 
 
            ierr = MPI_File_close(&fh);
            AssertThrowMPI(ierr);
          }
      } // if (mpisize > 1)
    else
#endif
      {
        (void)global_first_cell;
        (void)global_num_cells;
        (void)mpi_communicator;
 
        //
        // ---------- Fixed size data ----------
        //
        {
          const std::string fname_fixed =
            std::string(file_basename) + "_fixed.data";
 
          std::ofstream file(fname_fixed, std::ios::binary | std::ios::out);
          AssertThrow(file.fail() == false, ExcIO());
 
          // Write header data.
          file.write(reinterpret_cast<const char *>(
                       sizes_fixed_cumulative.data()),
                     sizes_fixed_cumulative.size() * sizeof(unsigned int));
 
          // Write packed data.
          file.write(reinterpret_cast<const char *>(src_data_fixed.data()),
                     src_data_fixed.size() * sizeof(char));
        }
 
        //
        // ---------- Variable size data ----------
        //
        if (variable_size_data_stored)
          {
            const std::string fname_variable =
              std::string(file_basename) + "_variable.data";
 
            std::ofstream file(fname_variable,
                               std::ios::binary | std::ios::out);
            AssertThrow(file.fail() == false, ExcIO());
 
            // Write header data.
            file.write(reinterpret_cast<const char *>(
                         src_sizes_variable.data()),
                       src_sizes_variable.size() * sizeof(int));
 
            // Write packed data.
            file.write(reinterpret_cast<const char *>(src_data_variable.data()),
                       src_data_variable.size() * sizeof(char));
          }
      }
  }
 
 
  template <int dim, int spacedim>
  DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
  void CellAttachedDataSerializer<dim, spacedim>::load(
    const unsigned int global_first_cell,
    const unsigned int global_num_cells,
    const unsigned int local_num_cells,
    const std::string &file_basename,
    const unsigned int n_attached_deserialize_fixed,
    const unsigned int n_attached_deserialize_variable,
    const MPI_Comm    &mpi_communicator)
  {
    Assert(dest_data_fixed.empty(),
           ExcMessage("Previously loaded data has not been released yet!"));
 
    variable_size_data_stored = (n_attached_deserialize_variable > 0);
 
#ifdef DEAL_II_WITH_MPI
    // Large fractions of this function have been copied from
    // DataOutInterface::write_vtu_in_parallel.
    // TODO: Write general MPIIO interface.
 
    const unsigned int mpisize =
      Utilities::MPI::n_mpi_processes(mpi_communicator);
 
    if (mpisize > 1)
      {
        //
        // ---------- Fixed size data ----------
        //
        {
          const std::string fname_fixed =
            std::string(file_basename) + "_fixed.data";
 
          MPI_Info info;
          int      ierr = MPI_Info_create(&info);
          AssertThrowMPI(ierr);
 
          MPI_File fh;
          ierr = MPI_File_open(
            mpi_communicator, fname_fixed.c_str(), MPI_MODE_RDONLY, info, &fh);
          AssertThrowMPI(ierr);
 
          ierr = MPI_Info_free(&info);
          AssertThrowMPI(ierr);
 
          // Read cumulative sizes from file.
          // Since all processors need the same information about the data
          // sizes, let each of them retrieve it by reading from the same
          // location in the file.
          sizes_fixed_cumulative.resize(1 + n_attached_deserialize_fixed +
                                        (variable_size_data_stored ? 1 : 0));
          ierr = Utilities::MPI::LargeCount::File_read_at_c(
            fh,
            0,
            sizes_fixed_cumulative.data(),
            sizes_fixed_cumulative.size(),
            MPI_UNSIGNED,
            MPI_STATUS_IGNORE);
          AssertThrowMPI(ierr);
 
          // Allocate sufficient memory.
          const unsigned int bytes_per_cell = sizes_fixed_cumulative.back();
          dest_data_fixed.resize(static_cast<size_t>(local_num_cells) *
                                 bytes_per_cell);
 
          // Read packed data from file simultaneously.
          const MPI_Offset size_header =
            sizes_fixed_cumulative.size() * sizeof(unsigned int);
 
          // Make sure we do the following computation in 64bit integers to be
          // able to handle 4GB+ files:
          const MPI_Offset my_global_file_position =
            size_header +
            static_cast<MPI_Offset>(global_first_cell) * bytes_per_cell;
 
          ierr =
            Utilities::MPI::LargeCount::File_read_at_c(fh,
                                                       my_global_file_position,
                                                       dest_data_fixed.data(),
                                                       dest_data_fixed.size(),
                                                       MPI_BYTE,
                                                       MPI_STATUS_IGNORE);
          AssertThrowMPI(ierr);
 
 
          ierr = MPI_File_close(&fh);
          AssertThrowMPI(ierr);
        }
 
        //
        // ---------- Variable size data ----------
        //
        if (variable_size_data_stored)
          {
            const std::string fname_variable =
              std::string(file_basename) + "_variable.data";
 
            MPI_Info info;
            int      ierr = MPI_Info_create(&info);
            AssertThrowMPI(ierr);
 
            MPI_File fh;
            ierr = MPI_File_open(mpi_communicator,
                                 fname_variable.c_str(),
                                 MPI_MODE_RDONLY,
                                 info,
                                 &fh);
            AssertThrowMPI(ierr);
 
            ierr = MPI_Info_free(&info);
            AssertThrowMPI(ierr);
 
            // Read sizes of all locally owned cells.
            dest_sizes_variable.resize(local_num_cells);
 
            const MPI_Offset my_global_file_position_sizes =
              static_cast<MPI_Offset>(global_first_cell) * sizeof(unsigned int);
 
            ierr = Utilities::MPI::LargeCount::File_read_at_c(
              fh,
              my_global_file_position_sizes,
              dest_sizes_variable.data(),
              dest_sizes_variable.size(),
              MPI_INT,
              MPI_STATUS_IGNORE);
            AssertThrowMPI(ierr);
 
 
            // Compute my data size in bytes and compute prefix sum. We do this
            // in 64 bit to avoid overflow for files larger than 4 GB:
            const std::uint64_t size_on_proc =
              std::accumulate(dest_sizes_variable.begin(),
                              dest_sizes_variable.end(),
                              0ULL);
 
            std::uint64_t prefix_sum = 0;
            ierr                     = MPI_Exscan(&size_on_proc,
                              &prefix_sum,
                              1,
                              MPI_UINT64_T,
                              MPI_SUM,
                              mpi_communicator);
            AssertThrowMPI(ierr);
 
            const MPI_Offset my_global_file_position =
              static_cast<MPI_Offset>(global_num_cells) * sizeof(unsigned int) +
              prefix_sum;
 
            dest_data_variable.resize(size_on_proc);
 
            ierr = Utilities::MPI::LargeCount::File_read_at_c(
              fh,
              my_global_file_position,
              dest_data_variable.data(),
              dest_data_variable.size(),
              MPI_BYTE,
              MPI_STATUS_IGNORE);
            AssertThrowMPI(ierr);
 
            ierr = MPI_File_close(&fh);
            AssertThrowMPI(ierr);
          }
      }
    else // if (mpisize > 1)
#endif
      {
        (void)global_first_cell;
        (void)global_num_cells;
        (void)mpi_communicator;
 
        //
        // ---------- Fixed size data ----------
        //
        {
          const std::string fname_fixed =
            std::string(file_basename) + "_fixed.data";
 
          std::ifstream file(fname_fixed, std::ios::binary | std::ios::in);
          AssertThrow(file.fail() == false, ExcIO());
 
          sizes_fixed_cumulative.resize(1 + n_attached_deserialize_fixed +
                                        (variable_size_data_stored ? 1 : 0));
 
          // Read header data.
          file.read(reinterpret_cast<char *>(sizes_fixed_cumulative.data()),
                    sizes_fixed_cumulative.size() * sizeof(unsigned int));
 
          const unsigned int bytes_per_cell = sizes_fixed_cumulative.back();
          dest_data_fixed.resize(static_cast<size_t>(local_num_cells) *
                                 bytes_per_cell);
 
          // Read packed data.
          file.read(reinterpret_cast<char *>(dest_data_fixed.data()),
                    dest_data_fixed.size() * sizeof(char));
        }
 
        //
        // ---------- Variable size data ----------
        //
        if (variable_size_data_stored)
          {
            const std::string fname_variable =
              std::string(file_basename) + "_variable.data";
 
            std::ifstream file(fname_variable, std::ios::binary | std::ios::in);
            AssertThrow(file.fail() == false, ExcIO());
 
            // Read header data.
            dest_sizes_variable.resize(local_num_cells);
            file.read(reinterpret_cast<char *>(dest_sizes_variable.data()),
                      dest_sizes_variable.size() * sizeof(int));
 
            // Read packed data.
            const std::uint64_t size =
              std::accumulate(dest_sizes_variable.begin(),
                              dest_sizes_variable.end(),
                              0ULL);
            dest_data_variable.resize(size);
            file.read(reinterpret_cast<char *>(dest_data_variable.data()),
                      dest_data_variable.size() * sizeof(char));
          }
      }
  }
 
 
  template <int dim, int spacedim>
  DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
  void CellAttachedDataSerializer<dim, spacedim>::clear()
  {
    variable_size_data_stored = false;
 
    // free information about data sizes
    sizes_fixed_cumulative.clear();
    sizes_fixed_cumulative.shrink_to_fit();
 
    // free fixed size transfer data
    src_data_fixed.clear();
    src_data_fixed.shrink_to_fit();
 
    dest_data_fixed.clear();
    dest_data_fixed.shrink_to_fit();
 
    // free variable size transfer data
    src_sizes_variable.clear();
    src_sizes_variable.shrink_to_fit();
 
    src_data_variable.clear();
    src_data_variable.shrink_to_fit();
 
    dest_sizes_variable.clear();
    dest_sizes_variable.shrink_to_fit();
 
    dest_data_variable.clear();
    dest_data_variable.shrink_to_fit();
  }
 
} // namespace internal
 
// anonymous namespace for internal helper functions
namespace
{
  // return whether the given cell is
  // patch_level_1, i.e. determine
  // whether either all or none of
  // its children are further
  // refined. this function can only
  // be called for non-active cells.
  template <int dim, int spacedim>
  bool
  cell_is_patch_level_1(
    const TriaIterator<::CellAccessor<dim, spacedim>> &cell)
  {
    Assert(cell->is_active() == false, ExcInternalError());
 
    unsigned int n_active_children = 0;
    for (unsigned int i = 0; i < cell->n_children(); ++i)
      if (cell->child(i)->is_active())
        ++n_active_children;
 
    return (n_active_children == 0) ||
           (n_active_children == cell->n_children());
  }
 
 
 
  // return, whether a given @p cell will be
  // coarsened, which is the case if all
  // children are active and have their coarsen
  // flag set. In case only part of the coarsen
  // flags are set, remove them.
  template <int dim, int spacedim>
  bool
  cell_will_be_coarsened(
    const TriaIterator<::CellAccessor<dim, spacedim>> &cell)
  {
    // only cells with children should be
    // considered for coarsening
 
    if (cell->has_children())
      {
        unsigned int       children_to_coarsen = 0;
        const unsigned int n_children          = cell->n_children();
 
        for (unsigned int c = 0; c < n_children; ++c)
          if (cell->child(c)->is_active() && cell->child(c)->coarsen_flag_set())
            ++children_to_coarsen;
        if (children_to_coarsen == n_children)
          return true;
        else
          for (unsigned int c = 0; c < n_children; ++c)
            if (cell->child(c)->is_active())
              cell->child(c)->clear_coarsen_flag();
      }
    // no children, so no coarsening
    // possible. however, no children also
    // means that this cell will be in the same
    // state as if it had children and was
    // coarsened. So, what should we return -
    // false or true?
    // make sure we do not have to do this at
    // all...
    Assert(cell->has_children(), ExcInternalError());
    // ... and then simply return false
    return false;
  }
 
 
  // return, whether the face @p face_no of the
  // given @p cell will be refined after the
  // current refinement step, considering
  // refine and coarsen flags and considering
  // only those refinemnts that will be caused
  // by the neighboring cell.
 
  // this function is used on both active cells
  // and cells with children. on cells with
  // children it also of interest to know 'how'
  // the face will be refined. thus there is an
  // additional third argument @p
  // expected_face_ref_case returning just
  // that. be aware, that this variable will
  // only contain useful information if this
  // function is called for an active cell.
  //
  // thus, this is an internal function, users
  // should call one of the two alternatives
  // following below.
  template <int dim, int spacedim>
  bool
  face_will_be_refined_by_neighbor_internal(
    const TriaIterator<::CellAccessor<dim, spacedim>> &cell,
    const unsigned int                                       face_no,
    RefinementCase<dim - 1> &expected_face_ref_case)
  {
    // first of all: set the default value for
    // expected_face_ref_case, which is no
    // refinement at all
    expected_face_ref_case = RefinementCase<dim - 1>::no_refinement;
 
    const typename Triangulation<dim, spacedim>::cell_iterator neighbor =
      cell->neighbor(face_no);
 
    // If we are at the boundary, there is no
    // neighbor which could refine the face
    if (neighbor.state() != IteratorState::valid)
      return false;
 
    if (neighbor->has_children())
      {
        // if the neighbor is refined, it may be
        // coarsened. if so, then it won't refine
        // the face, no matter what else happens
        if (cell_will_be_coarsened(neighbor))
          return false;
        else
          // if the neighbor is refined, then it
          // is also refined at our current
          // face. It will stay so without
          // coarsening, so return true in that
          // case.
          {
            expected_face_ref_case = cell->face(face_no)->refinement_case();
            return true;
          }
      }
 
    // now, the neighbor is not refined, but
    // perhaps it will be
    const RefinementCase<dim> nb_ref_flag = neighbor->refine_flag_set();
    if (nb_ref_flag != RefinementCase<dim>::no_refinement)
      {
        // now we need to know, which of the
        // neighbors faces points towards us
        const unsigned int neighbor_neighbor = cell->neighbor_face_no(face_no);
        // check, whether the cell will be
        // refined in a way that refines our
        // face
        const RefinementCase<dim - 1> face_ref_case =
          GeometryInfo<dim>::face_refinement_case(
            nb_ref_flag,
            neighbor_neighbor,
            neighbor->face_orientation(neighbor_neighbor),
            neighbor->face_flip(neighbor_neighbor),
            neighbor->face_rotation(neighbor_neighbor));
        if (face_ref_case != RefinementCase<dim - 1>::no_refinement)
          {
            const typename Triangulation<dim, spacedim>::face_iterator
                      neighbor_face   = neighbor->face(neighbor_neighbor);
            const int this_face_index = cell->face_index(face_no);
 
            // there are still two basic
            // possibilities here: the neighbor
            // might be coarser or as coarse
            // as we are
            if (neighbor_face->index() == this_face_index)
              // the neighbor is as coarse as
              // we are and will be refined at
              // the face of consideration, so
              // return true
              {
                expected_face_ref_case = face_ref_case;
                return true;
              }
            else
              {
                // the neighbor is coarser.
                // this is the most complicated
                // case. It might be, that the
                // neighbor's face will be
                // refined, but that we will
                // not see this, as we are
                // refined in a similar way.
 
                // so, the neighbor's face must
                // have children. check, if our
                // cell's face is one of these
                // (it could also be a
                // grand_child)
                for (unsigned int c = 0; c < neighbor_face->n_children(); ++c)
                  if (neighbor_face->child_index(c) == this_face_index)
                    {
                      // if the flagged refine
                      // case of the face is a
                      // subset or the same as
                      // the current refine case,
                      // then the face, as seen
                      // from our cell, won't be
                      // refined by the neighbor
                      if ((neighbor_face->refinement_case() | face_ref_case) ==
                          neighbor_face->refinement_case())
                        return false;
                      else
                        {
                          // if we are active, we
                          // must be an
                          // anisotropic child
                          // and the coming
                          // face_ref_case is
                          // isotropic. Thus,
                          // from our cell we
                          // will see exactly the
                          // opposite refine case
                          // that the face has
                          // now...
                          Assert(
                            face_ref_case ==
                              RefinementCase<dim - 1>::isotropic_refinement,
                            ExcInternalError());
                          expected_face_ref_case =
                            ~neighbor_face->refinement_case();
                          return true;
                        }
                    }
 
                // so, obviously we were not
                // one of the children, but a
                // grandchild. This is only
                // possible in 3d.
                Assert(dim == 3, ExcInternalError());
                // In that case, however, no
                // matter what the neighbor
                // does, it won't be finer
                // after the next refinement
                // step.
                return false;
              }
          } // if face will be refined
      }     // if neighbor is flagged for refinement
 
    // no cases left, so the neighbor will not
    // refine the face
    return false;
  }
 
  // version of above function for both active
  // and non-active cells
  template <int dim, int spacedim>
  bool
  face_will_be_refined_by_neighbor(
    const TriaIterator<::CellAccessor<dim, spacedim>> &cell,
    const unsigned int                                       face_no)
  {
    RefinementCase<dim - 1> dummy = RefinementCase<dim - 1>::no_refinement;
    return face_will_be_refined_by_neighbor_internal(cell, face_no, dummy);
  }
 
  // version of above function for active cells
  // only. Additionally returning the refine
  // case (to come) of the face under
  // consideration
  template <int dim, int spacedim>
  bool
  face_will_be_refined_by_neighbor(
    const TriaActiveIterator<::CellAccessor<dim, spacedim>> &cell,
    const unsigned int                                             face_no,
    RefinementCase<dim - 1> &expected_face_ref_case)
  {
    return face_will_be_refined_by_neighbor_internal(cell,
                                                     face_no,
                                                     expected_face_ref_case);
  }
 
 
 
  template <int dim, int spacedim>
  bool
  satisfies_level1_at_vertex_rule(
    const Triangulation<dim, spacedim> &triangulation)
  {
    std::vector<unsigned int> min_adjacent_cell_level(
      triangulation.n_vertices(), triangulation.n_levels());
    std::vector<unsigned int> max_adjacent_cell_level(
      triangulation.n_vertices(), 0);
 
    for (const auto &cell : triangulation.active_cell_iterators())
      for (const unsigned int v : cell->vertex_indices())
        {
          min_adjacent_cell_level[cell->vertex_index(v)] =
            std::min<unsigned int>(
              min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
          max_adjacent_cell_level[cell->vertex_index(v)] =
            std::max<unsigned int>(
              min_adjacent_cell_level[cell->vertex_index(v)], cell->level());
        }
 
    for (unsigned int k = 0; k < triangulation.n_vertices(); ++k)
      if (triangulation.vertex_used(k))
        if (max_adjacent_cell_level[k] - min_adjacent_cell_level[k] > 1)
          return false;
    return true;
  }
 
 
 
  template <int dim, int spacedim>
  unsigned int
  middle_vertex_index(
    const typename Triangulation<dim, spacedim>::line_iterator &line)
  {
    if (line->has_children())
      return line->child(0)->vertex_index(1);
    return numbers::invalid_unsigned_int;
  }
 
 
  template <int dim, int spacedim>
  unsigned int
  middle_vertex_index(
    const typename Triangulation<dim, spacedim>::quad_iterator &quad)
  {
    switch (static_cast<unsigned char>(quad->refinement_case()))
      {
        case RefinementCase<2>::cut_x:
          return middle_vertex_index<dim, spacedim>(quad->child(0)->line(1));
          break;
        case RefinementCase<2>::cut_y:
          return middle_vertex_index<dim, spacedim>(quad->child(0)->line(3));
          break;
        case RefinementCase<2>::cut_xy:
          return quad->child(0)->vertex_index(3);
          break;
        default:
          break;
      }
    return numbers::invalid_unsigned_int;
  }
 
 
  template <int dim, int spacedim>
  unsigned int
  middle_vertex_index(
    const typename Triangulation<dim, spacedim>::hex_iterator &hex)
  {
    switch (static_cast<unsigned char>(hex->refinement_case()))
      {
        case RefinementCase<3>::cut_x:
          return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(1));
          break;
        case RefinementCase<3>::cut_y:
          return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(3));
          break;
        case RefinementCase<3>::cut_z:
          return middle_vertex_index<dim, spacedim>(hex->child(0)->quad(5));
          break;
        case RefinementCase<3>::cut_xy:
          return middle_vertex_index<dim, spacedim>(hex->child(0)->line(11));
          break;
        case RefinementCase<3>::cut_xz:
          return middle_vertex_index<dim, spacedim>(hex->child(0)->line(5));
          break;
        case RefinementCase<3>::cut_yz:
          return middle_vertex_index<dim, spacedim>(hex->child(0)->line(7));
          break;
        case RefinementCase<3>::cut_xyz:
          return hex->child(0)->vertex_index(7);
          break;
        default:
          break;
      }
    return numbers::invalid_unsigned_int;
  }
 
 
  template <class TRIANGULATION>
  inline typename TRIANGULATION::DistortedCellList
  collect_distorted_coarse_cells(const TRIANGULATION &)
  {
    return typename TRIANGULATION::DistortedCellList();
  }
 
 
 
  template <int dim>
  inline typename Triangulation<dim, dim>::DistortedCellList
  collect_distorted_coarse_cells(const Triangulation<dim, dim> &triangulation)
  {
    typename Triangulation<dim, dim>::DistortedCellList distorted_cells;
    for (const auto &cell : triangulation.cell_iterators_on_level(0))
      {
        Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
        for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
          vertices[i] = cell->vertex(i);
 
        Tensor<0, dim> determinants[GeometryInfo<dim>::vertices_per_cell];
        GeometryInfo<dim>::alternating_form_at_vertices(vertices, determinants);
 
        for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
          if (determinants[i] <=
              1e-9 * Utilities::fixed_power<dim>(cell->diameter()))
            {
              distorted_cells.distorted_cells.push_back(cell);
              break;
            }
      }
 
    return distorted_cells;
  }
 
 
  template <int dim>
  bool
  has_distorted_children(
    const typename Triangulation<dim, dim>::cell_iterator &cell)
  {
    Assert(cell->has_children(), ExcInternalError());
 
    for (unsigned int c = 0; c < cell->n_children(); ++c)
      {
        Point<dim> vertices[GeometryInfo<dim>::vertices_per_cell];
        for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
          vertices[i] = cell->child(c)->vertex(i);
 
        Tensor<0, dim> determinants[GeometryInfo<dim>::vertices_per_cell];
        GeometryInfo<dim>::alternating_form_at_vertices(vertices, determinants);
 
        for (const unsigned int i : GeometryInfo<dim>::vertex_indices())
          if (determinants[i] <=
              1e-9 * Utilities::fixed_power<dim>(cell->child(c)->diameter()))
            return true;
      }
 
    return false;
  }
 
 
  template <int dim, int spacedim>
  bool
  has_distorted_children(
    const typename Triangulation<dim, spacedim>::cell_iterator &)
  {
    return false;
  }
 
 
  template <int dim, int spacedim>
  void
  update_periodic_face_map_recursively(
    const typename Triangulation<dim, spacedim>::cell_iterator &cell_1,
    const typename Triangulation<dim, spacedim>::cell_iterator &cell_2,
    unsigned int                                                n_face_1,
    unsigned int                                                n_face_2,
    const unsigned char                                         orientation,
    typename std::map<
      std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
                unsigned int>,
      std::pair<std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
                          unsigned int>,
                unsigned char>> &periodic_face_map)
  {
    using FaceIterator = typename Triangulation<dim, spacedim>::face_iterator;
    const FaceIterator face_1 = cell_1->face(n_face_1);
    const FaceIterator face_2 = cell_2->face(n_face_2);
 
    const unsigned char inverse_orientation =
      face_1->reference_cell().get_inverse_combined_orientation(orientation);
 
#ifdef DEBUG
    const auto [face_orientation, face_rotation, face_flip] =
      internal::split_face_orientation(orientation);
 
    Assert((dim != 1) || (face_orientation == true && face_flip == false &&
                          face_rotation == false),
           ExcMessage("The supplied orientation "
                      "(face_orientation, face_flip, face_rotation) "
                      "is invalid for 1d"));
 
    Assert((dim != 2) || (face_flip == false && face_rotation == false),
           ExcMessage("The supplied orientation "
                      "(face_orientation, face_flip, face_rotation) "
                      "is invalid for 2d"));
#endif
 
    Assert(face_1 != face_2, ExcMessage("face_1 and face_2 are equal!"));
 
    Assert(face_1->at_boundary() && face_2->at_boundary(),
           ExcMessage("Periodic faces must be on the boundary"));
 
    // Check if the requirement that each edge can only have at most one hanging
    // node, and as a consequence neighboring cells can differ by at most
    // one refinement level is enforced. In 1d, there are no hanging nodes and
    // so neighboring cells can differ by more than one refinement level.
    Assert(dim == 1 || std::abs(cell_1->level() - cell_2->level()) < 2,
           ExcInternalError());
 
    // insert periodic face pair for both cells
    using CellFace =
      std::pair<typename Triangulation<dim, spacedim>::cell_iterator,
                unsigned int>;
    const CellFace                           cell_face_1(cell_1, n_face_1);
    const CellFace                           cell_face_2(cell_2, n_face_2);
    const std::pair<CellFace, unsigned char> cell_face_orientation_2(
      cell_face_2, orientation);
 
    const std::pair<CellFace, std::pair<CellFace, unsigned char>>
      periodic_faces(cell_face_1, cell_face_orientation_2);
 
    // Only one periodic neighbor is allowed
    Assert(periodic_face_map.count(cell_face_1) == 0, ExcInternalError());
    periodic_face_map.insert(periodic_faces);
 
    if (dim == 1)
      {
        if (cell_1->has_children())
          {
            if (cell_2->has_children())
              {
                update_periodic_face_map_recursively<dim, spacedim>(
                  cell_1->child(n_face_1),
                  cell_2->child(n_face_2),
                  n_face_1,
                  n_face_2,
                  orientation,
                  periodic_face_map);
              }
            else // only face_1 has children
              {
                update_periodic_face_map_recursively<dim, spacedim>(
                  cell_1->child(n_face_1),
                  cell_2,
                  n_face_1,
                  n_face_2,
                  orientation,
                  periodic_face_map);
              }
          }
      }
    else // dim == 2 || dim == 3
      {
        if (cell_1->has_children())
          {
            if (cell_2->has_children())
              {
                // In the case that both faces have children, we loop over all
                // children and apply update_periodic_face_map_recursively
                // recursively:
 
                Assert(face_1->n_children() ==
                           GeometryInfo<dim>::max_children_per_face &&
                         face_2->n_children() ==
                           GeometryInfo<dim>::max_children_per_face,
                       ExcNotImplemented());
 
                const auto reference_cell = cell_1->reference_cell();
 
                for (unsigned int i = 0;
                     i < GeometryInfo<dim>::max_children_per_face;
                     ++i)
                  {
                    // Lookup the index for the second face
                    const unsigned int j =
                      reference_cell.standard_to_real_face_vertex(
                        i, n_face_1, inverse_orientation);
 
                    // find subcell ids that belong to the subface indices
                    unsigned int child_cell_1 =
                      GeometryInfo<dim>::child_cell_on_face(
                        cell_1->refinement_case(),
                        n_face_1,
                        i,
                        cell_1->face_orientation(n_face_1),
                        cell_1->face_flip(n_face_1),
                        cell_1->face_rotation(n_face_1),
                        face_1->refinement_case());
                    unsigned int child_cell_2 =
                      GeometryInfo<dim>::child_cell_on_face(
                        cell_2->refinement_case(),
                        n_face_2,
                        j,
                        cell_2->face_orientation(n_face_2),
                        cell_2->face_flip(n_face_2),
                        cell_2->face_rotation(n_face_2),
                        face_2->refinement_case());
 
                    Assert(cell_1->child(child_cell_1)->face(n_face_1) ==
                             face_1->child(i),
                           ExcInternalError());
                    Assert(cell_2->child(child_cell_2)->face(n_face_2) ==
                             face_2->child(j),
                           ExcInternalError());
 
                    // precondition: subcell has the same orientation as cell
                    // (so that the face numbers coincide) recursive call
                    update_periodic_face_map_recursively<dim, spacedim>(
                      cell_1->child(child_cell_1),
                      cell_2->child(child_cell_2),
                      n_face_1,
                      n_face_2,
                      orientation,
                      periodic_face_map);
                  }
              }
            else // only face_1 has children
              {
                for (unsigned int i = 0;
                     i < GeometryInfo<dim>::max_children_per_face;
                     ++i)
                  {
                    // find subcell ids that belong to the subface indices
                    unsigned int child_cell_1 =
                      GeometryInfo<dim>::child_cell_on_face(
                        cell_1->refinement_case(),
                        n_face_1,
                        i,
                        cell_1->face_orientation(n_face_1),
                        cell_1->face_flip(n_face_1),
                        cell_1->face_rotation(n_face_1),
                        face_1->refinement_case());
 
                    // recursive call
                    update_periodic_face_map_recursively<dim, spacedim>(
                      cell_1->child(child_cell_1),
                      cell_2,
                      n_face_1,
                      n_face_2,
                      orientation,
                      periodic_face_map);
                  }
              }
          }
      }
  }
 
 
} // end of anonymous namespace
 
 
namespace internal
{
  namespace TriangulationImplementation
  {
    // make sure that if in the following we
    // write Triangulation<dim,spacedim>
    // we mean the *class*
    // ::Triangulation, not the
    // enclosing namespace
    // internal::TriangulationImplementation
    using ::Triangulation;
 
    DeclException1(ExcInternalErrorOnCell,
                   int,
                   << "Something went wrong upon construction of cell "
                   << arg1);
    DeclException1(ExcCellHasNegativeMeasure,
                   int,
                   << "Cell " << arg1
                   << " has negative measure. This typically "
                   << "indicates some distortion in the cell, or a mistakenly "
                   << "swapped pair of vertices in the input to "
                   << "Triangulation::create_triangulation().");
    DeclException3(ExcInvalidVertexIndex,
                   int,
                   int,
                   int,
                   << "Error while creating cell " << arg1
                   << ": the vertex index " << arg2 << " must be between 0 and "
                   << arg3 << '.');
    DeclException3(
      ExcInteriorLineCantBeBoundary,
      int,
      int,
      types::boundary_id,
      << "The input data for creating a triangulation contained "
      << "information about a line with indices " << arg1 << " and " << arg2
      << " that is described to have boundary indicator "
      << static_cast<int>(arg3)
      << ". However, this is an internal line not located on the "
      << "boundary. You cannot assign a boundary indicator to it." << std::endl
      << std::endl
      << "If this happened at a place where you call "
      << "Triangulation::create_triangulation() yourself, you need "
      << "to check the SubCellData object you pass to this function."
      << std::endl
      << std::endl
      << "If this happened in a place where you are reading a mesh "
      << "from a file, then you need to investigate why such a line "
      << "ended up in the input file. A typical case is a geometry "
      << "that consisted of multiple parts and for which the mesh "
      << "generator program assumes that the interface between "
      << "two parts is a boundary when that isn't supposed to be "
      << "the case, or where the mesh generator simply assigns "
      << "'geometry indicators' to lines at the perimeter of "
      << "a part that are not supposed to be interpreted as "
      << "'boundary indicators'.");
    DeclException5(
      ExcInteriorQuadCantBeBoundary,
      int,
      int,
      int,
      int,
      types::boundary_id,
      << "The input data for creating a triangulation contained "
      << "information about a quad with indices " << arg1 << ", " << arg2
      << ", " << arg3 << ", and " << arg4
      << " that is described to have boundary indicator "
      << static_cast<int>(arg5)
      << ". However, this is an internal quad not located on the "
      << "boundary. You cannot assign a boundary indicator to it." << std::endl
      << std::endl
      << "If this happened at a place where you call "
      << "Triangulation::create_triangulation() yourself, you need "
      << "to check the SubCellData object you pass to this function."
      << std::endl
      << std::endl
      << "If this happened in a place where you are reading a mesh "
      << "from a file, then you need to investigate why such a quad "
      << "ended up in the input file. A typical case is a geometry "
      << "that consisted of multiple parts and for which the mesh "
      << "generator program assumes that the interface between "
      << "two parts is a boundary when that isn't supposed to be "
      << "the case, or where the mesh generator simply assigns "
      << "'geometry indicators' to quads at the surface of "
      << "a part that are not supposed to be interpreted as "
      << "'boundary indicators'.");
    DeclException2(
      ExcMultiplySetLineInfoOfLine,
      int,
      int,
      << "In SubCellData the line info of the line with vertex indices " << arg1
      << " and " << arg2 << " appears more than once. "
      << "This is not allowed.");
    DeclException3(
      ExcInconsistentLineInfoOfLine,
      int,
      int,
      std::string,
      << "In SubCellData the line info of the line with vertex indices " << arg1
      << " and " << arg2 << " appears multiple times with different (valid) "
      << arg3 << ". This is not allowed.");
    DeclException5(
      ExcInconsistentQuadInfoOfQuad,
      int,
      int,
      int,
      int,
      std::string,
      << "In SubCellData the quad info of the quad with line indices " << arg1
      << ", " << arg2 << ", " << arg3 << " and " << arg4
      << " appears multiple times with different (valid) " << arg5
      << ". This is not allowed.");
 
    /*
     * Reserve space for TriaFaces. Details:
     *
     * Reserve space for line_orientations.
     *
     * @note Used only for dim=3.
     */
    void
    reserve_space(TriaFaces         &tria_faces,
                  const unsigned int new_quads_in_pairs,
                  const unsigned int new_quads_single)
    {
      AssertDimension(tria_faces.dim, 3);
 
      Assert(new_quads_in_pairs % 2 == 0, ExcInternalError());
 
      unsigned int next_free_single = 0;
      unsigned int next_free_pair   = 0;
 
      // count the number of objects, of unused single objects and of
      // unused pairs of objects
      unsigned int n_quads          = 0;
      unsigned int n_unused_pairs   = 0;
      unsigned int n_unused_singles = 0;
      for (unsigned int i = 0; i < tria_faces.quads.used.size(); ++i)
        {
          if (tria_faces.quads.used[i])
            ++n_quads;
          else if (i + 1 < tria_faces.quads.used.size())
            {
              if (tria_faces.quads.used[i + 1])
                {
                  ++n_unused_singles;
                  if (next_free_single == 0)
                    next_free_single = i;
                }
              else
                {
                  ++n_unused_pairs;
                  if (next_free_pair == 0)
                    next_free_pair = i;
                  ++i;
                }
            }
          else
            ++n_unused_singles;
        }
      Assert(n_quads + 2 * n_unused_pairs + n_unused_singles ==
               tria_faces.quads.used.size(),
             ExcInternalError());
      (void)n_quads;
 
      // how many single quads are needed in addition to n_unused_quads?
      const int additional_single_quads = new_quads_single - n_unused_singles;
 
      unsigned int new_size =
        tria_faces.quads.used.size() + new_quads_in_pairs - 2 * n_unused_pairs;
      if (additional_single_quads > 0)
        new_size += additional_single_quads;
 
      // see above...
      if (new_size > tria_faces.quads.n_objects())
        {
          // reserve the field of the derived class
          tria_faces.quads_line_orientations.resize(
            new_size * GeometryInfo<3>::lines_per_face, true);
 
          auto &q_is_q = tria_faces.quad_is_quadrilateral;
          q_is_q.reserve(new_size);
          q_is_q.insert(q_is_q.end(), new_size - q_is_q.size(), true);
        }
    }
 
 
 
    void
    reserve_space(TriaLevel         &tria_level,
                  const unsigned int total_cells,
                  const unsigned int dimension,
                  const unsigned int space_dimension,
                  const bool         tetraheder_in_mesh = false)
    {
      // we need space for total_cells cells. Maybe we have more already
      // with those cells which are unused, so only allocate new space if
      // needed.
      //
      // note that all arrays should have equal sizes (checked by
      // @p{monitor_memory}
      if (total_cells > tria_level.refine_flags.size())
        {
          tria_level.refine_flags.reserve(total_cells);
          tria_level.refine_flags.insert(tria_level.refine_flags.end(),
                                         total_cells -
                                           tria_level.refine_flags.size(),
                                         /*RefinementCase::no_refinement=*/0);
 
          if (tetraheder_in_mesh)
            {
              tria_level.refine_choice.reserve(total_cells);
              tria_level.refine_choice.insert(
                tria_level.refine_choice.end(),
                total_cells - tria_level.refine_choice.size(),
                static_cast<char>(
                  IsotropicRefinementChoice::isotropic_refinement));
            }
 
          tria_level.coarsen_flags.reserve(total_cells);
          tria_level.coarsen_flags.insert(tria_level.coarsen_flags.end(),
                                          total_cells -
                                            tria_level.coarsen_flags.size(),
                                          false);
 
          tria_level.active_cell_indices.reserve(total_cells);
          tria_level.active_cell_indices.insert(
            tria_level.active_cell_indices.end(),
            total_cells - tria_level.active_cell_indices.size(),
            numbers::invalid_unsigned_int);
 
          tria_level.subdomain_ids.reserve(total_cells);
          tria_level.subdomain_ids.insert(tria_level.subdomain_ids.end(),
                                          total_cells -
                                            tria_level.subdomain_ids.size(),
                                          0);
 
          tria_level.level_subdomain_ids.reserve(total_cells);
          tria_level.level_subdomain_ids.insert(
            tria_level.level_subdomain_ids.end(),
            total_cells - tria_level.level_subdomain_ids.size(),
            0);
 
          tria_level.global_active_cell_indices.reserve(total_cells);
          tria_level.global_active_cell_indices.insert(
            tria_level.global_active_cell_indices.end(),
            total_cells - tria_level.global_active_cell_indices.size(),
            numbers::invalid_dof_index);
 
          tria_level.global_level_cell_indices.reserve(total_cells);
          tria_level.global_level_cell_indices.insert(
            tria_level.global_level_cell_indices.end(),
            total_cells - tria_level.global_level_cell_indices.size(),
            numbers::invalid_dof_index);
 
          if (dimension == space_dimension - 1)
            {
              tria_level.direction_flags.reserve(total_cells);
              tria_level.direction_flags.insert(
                tria_level.direction_flags.end(),
                total_cells - tria_level.direction_flags.size(),
                true);
            }
          else
            tria_level.direction_flags.clear();
 
          tria_level.parents.reserve((total_cells + 1) / 2);
          tria_level.parents.insert(tria_level.parents.end(),
                                    (total_cells + 1) / 2 -
                                      tria_level.parents.size(),
                                    -1);
 
          tria_level.neighbors.reserve(total_cells * (2 * dimension));
          tria_level.neighbors.insert(tria_level.neighbors.end(),
                                      total_cells * (2 * dimension) -
                                        tria_level.neighbors.size(),
                                      std::make_pair(-1, -1));
 
          if (tria_level.dim == 2 || tria_level.dim == 3)
            {
              const unsigned int max_faces_per_cell = 2 * dimension;
              tria_level.face_orientations.resize(total_cells *
                                                  max_faces_per_cell);
 
              tria_level.reference_cell.reserve(total_cells);
              tria_level.reference_cell.insert(
                tria_level.reference_cell.end(),
                total_cells - tria_level.reference_cell.size(),
                tria_level.dim == 2 ? ReferenceCells::Quadrilateral :
                                      ReferenceCells::Hexahedron);
            }
        }
    }
 
 
 
    DeclException2(ExcMemoryInexact,
                   int,
                   int,
                   << "The containers have sizes " << arg1 << " and " << arg2
                   << ", which is not as expected.");
 
    void
    monitor_memory(const TriaLevel   &tria_level,
                   const unsigned int true_dimension)
    {
      (void)tria_level;
      (void)true_dimension;
      Assert(2 * true_dimension * tria_level.refine_flags.size() ==
               tria_level.neighbors.size(),
             ExcMemoryInexact(tria_level.refine_flags.size(),
                              tria_level.neighbors.size()));
      Assert(2 * true_dimension * tria_level.coarsen_flags.size() ==
               tria_level.neighbors.size(),
             ExcMemoryInexact(tria_level.coarsen_flags.size(),
                              tria_level.neighbors.size()));
    }
 
 
 
    void
    reserve_space(TriaObjects       &tria_objects,
                  const unsigned int new_objects_in_pairs,
                  const unsigned int new_objects_single = 0)
    {
      if (tria_objects.structdim <= 2)
        {
          Assert(new_objects_in_pairs % 2 == 0, ExcInternalError());
 
          tria_objects.next_free_single               = 0;
          tria_objects.next_free_pair                 = 0;
          tria_objects.reverse_order_next_free_single = false;
 
          // count the number of objects, of unused single objects and of
          // unused pairs of objects
          unsigned int n_objects        = 0;
          unsigned int n_unused_pairs   = 0;
          unsigned int n_unused_singles = 0;
          for (unsigned int i = 0; i < tria_objects.used.size(); ++i)
            {
              if (tria_objects.used[i])
                ++n_objects;
              else if (i + 1 < tria_objects.used.size())
                {
                  if (tria_objects.used[i + 1])
                    {
                      ++n_unused_singles;
                      if (tria_objects.next_free_single == 0)
                        tria_objects.next_free_single = i;
                    }
                  else
                    {
                      ++n_unused_pairs;
                      if (tria_objects.next_free_pair == 0)
                        tria_objects.next_free_pair = i;
                      ++i;
                    }
                }
              else
                ++n_unused_singles;
            }
          Assert(n_objects + 2 * n_unused_pairs + n_unused_singles ==
                   tria_objects.used.size(),
                 ExcInternalError());
          (void)n_objects;
 
          // how many single objects are needed in addition to
          // n_unused_objects?
          const int additional_single_objects =
            new_objects_single - n_unused_singles;
 
          unsigned int new_size = tria_objects.used.size() +
                                  new_objects_in_pairs - 2 * n_unused_pairs;
          if (additional_single_objects > 0)
            new_size += additional_single_objects;
 
          // only allocate space if necessary
          if (new_size > tria_objects.n_objects())
            {
              const unsigned int max_faces_per_cell =
                2 * tria_objects.structdim;
              const unsigned int max_children_per_cell =
                1 << tria_objects.structdim;
 
              tria_objects.cells.reserve(new_size * max_faces_per_cell);
              tria_objects.cells.insert(tria_objects.cells.end(),
                                        (new_size - tria_objects.n_objects()) *
                                          max_faces_per_cell,
                                        -1);
 
              tria_objects.used.reserve(new_size);
              tria_objects.used.insert(tria_objects.used.end(),
                                       new_size - tria_objects.used.size(),
                                       false);
 
              tria_objects.user_flags.reserve(new_size);
              tria_objects.user_flags.insert(tria_objects.user_flags.end(),
                                             new_size -
                                               tria_objects.user_flags.size(),
                                             false);
 
              const unsigned int factor = max_children_per_cell / 2;
              tria_objects.children.reserve(factor * new_size);
              tria_objects.children.insert(tria_objects.children.end(),
                                           factor * new_size -
                                             tria_objects.children.size(),
                                           -1);
 
              if (tria_objects.structdim > 1)
                {
                  tria_objects.refinement_cases.reserve(new_size);
                  tria_objects.refinement_cases.insert(
                    tria_objects.refinement_cases.end(),
                    new_size - tria_objects.refinement_cases.size(),
                    /*RefinementCase::no_refinement=*/0);
                }
 
              // first reserve, then resize. Otherwise the std library can
              // decide to allocate more entries.
              tria_objects.boundary_or_material_id.reserve(new_size);
              tria_objects.boundary_or_material_id.resize(new_size);
 
              tria_objects.user_data.reserve(new_size);
              tria_objects.user_data.resize(new_size);
 
              tria_objects.manifold_id.reserve(new_size);
              tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
                                              new_size -
                                                tria_objects.manifold_id.size(),
                                              numbers::flat_manifold_id);
            }
 
          if (n_unused_singles == 0)
            {
              tria_objects.next_free_single               = new_size - 1;
              tria_objects.reverse_order_next_free_single = true;
            }
        }
      else
        {
          const unsigned int new_hexes = new_objects_in_pairs;
 
          const unsigned int new_size =
            new_hexes + std::count(tria_objects.used.begin(),
                                   tria_objects.used.end(),
                                   true);
 
          // see above...
          if (new_size > tria_objects.n_objects())
            {
              const unsigned int max_faces_per_cell =
                2 * tria_objects.structdim;
 
              tria_objects.cells.reserve(new_size * max_faces_per_cell);
              tria_objects.cells.insert(tria_objects.cells.end(),
                                        (new_size - tria_objects.n_objects()) *
                                          max_faces_per_cell,
                                        -1);
 
              tria_objects.used.reserve(new_size);
              tria_objects.used.insert(tria_objects.used.end(),
                                       new_size - tria_objects.used.size(),
                                       false);
 
              tria_objects.user_flags.reserve(new_size);
              tria_objects.user_flags.insert(tria_objects.user_flags.end(),
                                             new_size -
                                               tria_objects.user_flags.size(),
                                             false);
 
              tria_objects.children.reserve(4 * new_size);
              tria_objects.children.insert(tria_objects.children.end(),
                                           4 * new_size -
                                             tria_objects.children.size(),
                                           -1);
 
              // for the following fields, we know exactly how many elements
              // we need, so first reserve then resize (resize itself, at least
              // with some compiler libraries, appears to round up the size it
              // actually reserves)
              tria_objects.boundary_or_material_id.reserve(new_size);
              tria_objects.boundary_or_material_id.resize(new_size);
 
              tria_objects.manifold_id.reserve(new_size);
              tria_objects.manifold_id.insert(tria_objects.manifold_id.end(),
                                              new_size -
                                                tria_objects.manifold_id.size(),
                                              numbers::flat_manifold_id);
 
              tria_objects.user_data.reserve(new_size);
              tria_objects.user_data.resize(new_size);
 
              tria_objects.refinement_cases.reserve(new_size);
              tria_objects.refinement_cases.insert(
                tria_objects.refinement_cases.end(),
                new_size - tria_objects.refinement_cases.size(),
                /*RefinementCase::no_refinement=*/0);
            }
          tria_objects.next_free_single = tria_objects.next_free_pair = 0;
        }
    }
 
 
 
    void
    monitor_memory(const TriaObjects &tria_object, const unsigned int)
    {
      Assert(tria_object.n_objects() == tria_object.used.size(),
             ExcMemoryInexact(tria_object.n_objects(),
                              tria_object.used.size()));
      Assert(tria_object.n_objects() == tria_object.user_flags.size(),
             ExcMemoryInexact(tria_object.n_objects(),
                              tria_object.user_flags.size()));
      Assert(tria_object.n_objects() ==
               tria_object.boundary_or_material_id.size(),
             ExcMemoryInexact(tria_object.n_objects(),
                              tria_object.boundary_or_material_id.size()));
      Assert(tria_object.n_objects() == tria_object.manifold_id.size(),
             ExcMemoryInexact(tria_object.n_objects(),
                              tria_object.manifold_id.size()));
      Assert(tria_object.n_objects() == tria_object.user_data.size(),
             ExcMemoryInexact(tria_object.n_objects(),
                              tria_object.user_data.size()));
 
      if (tria_object.structdim == 1)
        {
          Assert(1 * tria_object.n_objects() == tria_object.children.size(),
                 ExcMemoryInexact(tria_object.n_objects(),
                                  tria_object.children.size()));
        }
      else if (tria_object.structdim == 2)
        {
          Assert(2 * tria_object.n_objects() == tria_object.children.size(),
                 ExcMemoryInexact(tria_object.n_objects(),
                                  tria_object.children.size()));
        }
      else if (tria_object.structdim == 3)
        {
          Assert(4 * tria_object.n_objects() == tria_object.children.size(),
                 ExcMemoryInexact(tria_object.n_objects(),
                                  tria_object.children.size()));
        }
    }
 
 
 
    template <int dim, int spacedim>
    class Policy
    {
    public:
      virtual ~Policy() = default;
 
      virtual void
      update_neighbors(Triangulation<dim, spacedim> &tria) = 0;
 
      virtual void
      delete_children(
        Triangulation<dim, spacedim>                         &triangulation,
        typename Triangulation<dim, spacedim>::cell_iterator &cell,
        std::vector<unsigned int>                            &line_cell_count,
        std::vector<unsigned int> &quad_cell_count) = 0;
 
      virtual typename Triangulation<dim, spacedim>::DistortedCellList
      execute_refinement(Triangulation<dim, spacedim> &triangulation,
                         const bool check_for_distorted_cells) = 0;
 
      virtual void
      prevent_distorted_boundary_cells(
        Triangulation<dim, spacedim> &triangulation) = 0;
 
      virtual void
      prepare_refinement_dim_dependent(
        Triangulation<dim, spacedim> &triangulation) = 0;
 
      virtual bool
      coarsening_allowed(
        const typename Triangulation<dim, spacedim>::cell_iterator &cell) = 0;
 
      virtual std::unique_ptr<Policy<dim, spacedim>>
      clone() = 0;
    };
 
 
 
    template <int dim, int spacedim, typename T>
    class PolicyWrapper : public Policy<dim, spacedim>
    {
    public:
      void
      update_neighbors(Triangulation<dim, spacedim> &tria) override
      {
        T::update_neighbors(tria);
      }
 
      void
      delete_children(
        Triangulation<dim, spacedim>                         &tria,
        typename Triangulation<dim, spacedim>::cell_iterator &cell,
        std::vector<unsigned int>                            &line_cell_count,
        std::vector<unsigned int> &quad_cell_count) override
      {
        T::delete_children(tria, cell, line_cell_count, quad_cell_count);
      }
 
      typename Triangulation<dim, spacedim>::DistortedCellList
      execute_refinement(Triangulation<dim, spacedim> &triangulation,
                         const bool check_for_distorted_cells) override
      {
        return T::execute_refinement(triangulation, check_for_distorted_cells);
      }
 
      void
      prevent_distorted_boundary_cells(
        Triangulation<dim, spacedim> &triangulation) override
      {
        T::prevent_distorted_boundary_cells(triangulation);
      }
 
      void
      prepare_refinement_dim_dependent(
        Triangulation<dim, spacedim> &triangulation) override
      {
        T::prepare_refinement_dim_dependent(triangulation);
      }
 
      bool
      coarsening_allowed(
        const typename Triangulation<dim, spacedim>::cell_iterator &cell)
        override
      {
        return T::template coarsening_allowed<dim, spacedim>(cell);
      }
 
      std::unique_ptr<Policy<dim, spacedim>>
      clone() override
      {
        return std::make_unique<PolicyWrapper<dim, spacedim, T>>();
      }
    };
 
 
 
    struct Implementation
    {
      template <int dim, int spacedim>
      static void
      compute_number_cache_dim(
        const Triangulation<dim, spacedim>                    &triangulation,
        const unsigned int                                     level_objects,
        internal::TriangulationImplementation::NumberCache<1> &number_cache)
      {
        using line_iterator =
          typename Triangulation<dim, spacedim>::line_iterator;
 
        number_cache.n_levels = 0;
        if (level_objects > 0)
          // find the last level on which there are used cells
          for (unsigned int level = 0; level < level_objects; ++level)
            if (triangulation.begin(level) != triangulation.end(level))
              number_cache.n_levels = level + 1;
 
        // no cells at all?
        Assert(number_cache.n_levels > 0, ExcInternalError());
 
        //---------------------------------
        // update the number of lines on the different levels in the
        // cache
        number_cache.n_lines        = 0;
        number_cache.n_active_lines = 0;
 
        // for 1d, lines have levels so take count the objects per
        // level and globally
        if (dim == 1)
          {
            number_cache.n_lines_level.resize(number_cache.n_levels);
            number_cache.n_active_lines_level.resize(number_cache.n_levels);
 
            for (unsigned int level = 0; level < number_cache.n_levels; ++level)
              {
                // count lines on this level
                number_cache.n_lines_level[level]        = 0;
                number_cache.n_active_lines_level[level] = 0;
 
                line_iterator line = triangulation.begin_line(level),
                              endc =
                                (level == number_cache.n_levels - 1 ?
                                   line_iterator(triangulation.end_line()) :
                                   triangulation.begin_line(level + 1));
                for (; line != endc; ++line)
                  {
                    ++number_cache.n_lines_level[level];
                    if (line->has_children() == false)
                      ++number_cache.n_active_lines_level[level];
                  }
 
                // update total number of lines
                number_cache.n_lines += number_cache.n_lines_level[level];
                number_cache.n_active_lines +=
                  number_cache.n_active_lines_level[level];
              }
          }
        else
          {
            // for dim>1, there are no levels for lines
            number_cache.n_lines_level.clear();
            number_cache.n_active_lines_level.clear();
 
            line_iterator line = triangulation.begin_line(),
                          endc = triangulation.end_line();
            for (; line != endc; ++line)
              {
                ++number_cache.n_lines;
                if (line->has_children() == false)
                  ++number_cache.n_active_lines;
              }
          }
      }
 
      template <int dim, int spacedim>
      static void
      compute_number_cache_dim(
        const Triangulation<dim, spacedim>                    &triangulation,
        const unsigned int                                     level_objects,
        internal::TriangulationImplementation::NumberCache<2> &number_cache)
      {
        // update lines and n_levels in number_cache. since we don't
        // access any of these numbers, we can do this in the
        // background
        Threads::Task<void> update_lines = Threads::new_task(
          static_cast<
            void (*)(const Triangulation<dim, spacedim> &,
                     const unsigned int,
                     internal::TriangulationImplementation::NumberCache<1> &)>(
            &compute_number_cache_dim<dim, spacedim>),
          triangulation,
          level_objects,
          static_cast<internal::TriangulationImplementation::NumberCache<1> &>(
            number_cache));
 
        using quad_iterator =
          typename Triangulation<dim, spacedim>::quad_iterator;
 
        //---------------------------------
        // update the number of quads on the different levels in the
        // cache
        number_cache.n_quads        = 0;
        number_cache.n_active_quads = 0;
 
        // for 2d, quads have levels so take count the objects per
        // level and globally
        if (dim == 2)
          {
            // count the number of levels; the function we called above
            // on a separate Task for lines also does this and puts it into
            // number_cache.n_levels, but this datum may not yet be
            // available as we call the function on a separate task
            unsigned int n_levels = 0;
            if (level_objects > 0)
              // find the last level on which there are used cells
              for (unsigned int level = 0; level < level_objects; ++level)
                if (triangulation.begin(level) != triangulation.end(level))
                  n_levels = level + 1;
 
            number_cache.n_quads_level.resize(n_levels);
            number_cache.n_active_quads_level.resize(n_levels);
 
            for (unsigned int level = 0; level < n_levels; ++level)
              {
                // count quads on this level
                number_cache.n_quads_level[level]        = 0;
                number_cache.n_active_quads_level[level] = 0;
 
                quad_iterator quad = triangulation.begin_quad(level),
                              endc =
                                (level == n_levels - 1 ?
                                   quad_iterator(triangulation.end_quad()) :
                                   triangulation.begin_quad(level + 1));
                for (; quad != endc; ++quad)
                  {
                    ++number_cache.n_quads_level[level];
                    if (quad->has_children() == false)
                      ++number_cache.n_active_quads_level[level];
                  }
 
                // update total number of quads
                number_cache.n_quads += number_cache.n_quads_level[level];
                number_cache.n_active_quads +=
                  number_cache.n_active_quads_level[level];
              }
          }
        else
          {
            // for dim>2, there are no levels for quads
            number_cache.n_quads_level.clear();
            number_cache.n_active_quads_level.clear();
 
            quad_iterator quad = triangulation.begin_quad(),
                          endc = triangulation.end_quad();
            for (; quad != endc; ++quad)
              {
                ++number_cache.n_quads;
                if (quad->has_children() == false)
                  ++number_cache.n_active_quads;
              }
          }
 
        // wait for the background computation for lines
        update_lines.join();
      }
 
      template <int dim, int spacedim>
      static void
      compute_number_cache_dim(
        const Triangulation<dim, spacedim>                    &triangulation,
        const unsigned int                                     level_objects,
        internal::TriangulationImplementation::NumberCache<3> &number_cache)
      {
        // update quads, lines and n_levels in number_cache. since we
        // don't access any of these numbers, we can do this in the
        // background
        Threads::Task<void> update_quads_and_lines = Threads::new_task(
          static_cast<
            void (*)(const Triangulation<dim, spacedim> &,
                     const unsigned int,
                     internal::TriangulationImplementation::NumberCache<2> &)>(
            &compute_number_cache_dim<dim, spacedim>),
          triangulation,
          level_objects,
          static_cast<internal::TriangulationImplementation::NumberCache<2> &>(
            number_cache));
 
        using hex_iterator =
          typename Triangulation<dim, spacedim>::hex_iterator;
 
        //---------------------------------
        // update the number of hexes on the different levels in the
        // cache
        number_cache.n_hexes        = 0;
        number_cache.n_active_hexes = 0;
 
        // for 3d, hexes have levels so take count the objects per
        // level and globally
        if (dim == 3)
          {
            // count the number of levels; the function we called
            // above on a separate Task for quads (recursively, via
            // the lines function) also does this and puts it into
            // number_cache.n_levels, but this datum may not yet be
            // available as we call the function on a separate task
            unsigned int n_levels = 0;
            if (level_objects > 0)
              // find the last level on which there are used cells
              for (unsigned int level = 0; level < level_objects; ++level)
                if (triangulation.begin(level) != triangulation.end(level))
                  n_levels = level + 1;
 
            number_cache.n_hexes_level.resize(n_levels);
            number_cache.n_active_hexes_level.resize(n_levels);
 
            for (unsigned int level = 0; level < n_levels; ++level)
              {
                // count hexes on this level
                number_cache.n_hexes_level[level]        = 0;
                number_cache.n_active_hexes_level[level] = 0;
 
                hex_iterator hex  = triangulation.begin_hex(level),
                             endc = (level == n_levels - 1 ?
                                       hex_iterator(triangulation.end_hex()) :
                                       triangulation.begin_hex(level + 1));
                for (; hex != endc; ++hex)
                  {
                    ++number_cache.n_hexes_level[level];
                    if (hex->has_children() == false)
                      ++number_cache.n_active_hexes_level[level];
                  }
 
                // update total number of hexes
                number_cache.n_hexes += number_cache.n_hexes_level[level];
                number_cache.n_active_hexes +=
                  number_cache.n_active_hexes_level[level];
              }
          }
        else
          {
            // for dim>3, there are no levels for hexes
            number_cache.n_hexes_level.clear();
            number_cache.n_active_hexes_level.clear();
 
            hex_iterator hex  = triangulation.begin_hex(),
                         endc = triangulation.end_hex();
            for (; hex != endc; ++hex)
              {
                ++number_cache.n_hexes;
                if (hex->has_children() == false)
                  ++number_cache.n_active_hexes;
              }
          }
 
        // wait for the background computation for quads
        update_quads_and_lines.join();
      }
 
 
      template <int dim, int spacedim>
      static void
      compute_number_cache(
        const Triangulation<dim, spacedim>                      &triangulation,
        const unsigned int                                       level_objects,
        internal::TriangulationImplementation::NumberCache<dim> &number_cache)
      {
        compute_number_cache_dim(triangulation, level_objects, number_cache);
 
        number_cache.active_cell_index_partitioner =
          std::make_shared<const Utilities::MPI::Partitioner>(
            triangulation.n_active_cells());
 
        number_cache.level_cell_index_partitioners.resize(
          triangulation.n_levels());
        for (unsigned int level = 0; level < triangulation.n_levels(); ++level)
          number_cache.level_cell_index_partitioners[level] =
            std::make_shared<const Utilities::MPI::Partitioner>(
              triangulation.n_cells(level));
      }
 
 
      template <int spacedim>
      static void
      update_neighbors(Triangulation<1, spacedim> &)
      {}
 
 
      template <int dim, int spacedim>
      static void
      update_neighbors(Triangulation<dim, spacedim> &triangulation)
      {
        // each face can be neighbored on two sides
        // by cells. according to the face's
        // intrinsic normal we define the left
        // neighbor as the one for which the face
        // normal points outward, and store that
        // one first; the second one is then
        // the right neighbor for which the
        // face normal points inward. This
        // information depends on the type of cell
        // and local number of face for the
        // 'standard ordering and orientation' of
        // faces and then on the face_orientation
        // information for the real mesh. Set up a
        // table to have fast access to those
        // offsets (0 for left and 1 for
        // right). Some of the values are invalid
        // as they reference too large face
        // numbers, but we just leave them at a
        // zero value.
        //
        // Note, that in 2d for lines as faces the
        // normal direction given in the
        // GeometryInfo class is not consistent. We
        // thus define here that the normal for a
        // line points to the right if the line
        // points upwards.
        //
        // There is one more point to
        // consider, however: if we have
        // dim<spacedim, then we may have
        // cases where cells are
        // inverted. In effect, both
        // cells think they are the left
        // neighbor of an edge, for
        // example, which leads us to
        // forget neighborship
        // information (a case that shows
        // this is
        // codim_one/hanging_nodes_02). We
        // store whether a cell is
        // inverted using the
        // direction_flag, so if a cell
        // has a false direction_flag,
        // then we need to invert our
        // selection whether we are a
        // left or right neighbor in all
        // following computations.
        //
        // first index:  dimension (minus 2)
        // second index: local face index
        // third index:  face_orientation (false and true)
        static const unsigned int left_right_offset[2][6][2] = {
          // quadrilateral
          {{0, 1},  // face 0, face_orientation = false and true
           {1, 0},  // face 1, face_orientation = false and true
           {1, 0},  // face 2, face_orientation = false and true
           {0, 1},  // face 3, face_orientation = false and true
           {0, 0},  // face 4, invalid face
           {0, 0}}, // face 5, invalid face
                    // hexahedron
          {{0, 1}, {1, 0}, {0, 1}, {1, 0}, {0, 1}, {1, 0}}};
 
        // now create a vector of the two active
        // neighbors (left and right) for each face
        // and fill it by looping over all cells. For
        // cases with anisotropic refinement and more
        // then one cell neighboring at a given side
        // of the face we will automatically get the
        // active one on the highest level as we loop
        // over cells from lower levels first.
        const typename Triangulation<dim, spacedim>::cell_iterator dummy;
        std::vector<typename Triangulation<dim, spacedim>::cell_iterator>
          adjacent_cells(2 * triangulation.n_raw_faces(), dummy);
 
        for (const auto &cell : triangulation.cell_iterators())
          for (auto f : cell->face_indices())
            {
              const typename Triangulation<dim, spacedim>::face_iterator face =
                cell->face(f);
 
              const unsigned int offset =
                (cell->direction_flag() ?
                   left_right_offset[dim - 2][f][cell->face_orientation(f)] :
                   1 -
                     left_right_offset[dim - 2][f][cell->face_orientation(f)]);
 
              adjacent_cells[2 * face->index() + offset] = cell;
 
              // if this cell is not refined, but the
              // face is, then we'll have to set our
              // cell as neighbor for the child faces
              // as well. Fortunately the normal
              // orientation of children will be just
              // the same.
              if (dim == 2)
                {
                  if (cell->is_active() && face->has_children())
                    {
                      adjacent_cells[2 * face->child(0)->index() + offset] =
                        cell;
                      adjacent_cells[2 * face->child(1)->index() + offset] =
                        cell;
                    }
                }
              else // -> dim == 3
                {
                  // We need the same as in 2d
                  // here. Furthermore, if the face is
                  // refined with cut_x or cut_y then
                  // those children again in the other
                  // direction, and if this cell is
                  // refined isotropically (along the
                  // face) then the neighbor will
                  // (probably) be refined as cut_x or
                  // cut_y along the face. For those
                  // neighboring children cells, their
                  // neighbor will be the current,
                  // inactive cell, as our children are
                  // too fine to be neighbors. Catch that
                  // case by also acting on inactive
                  // cells with isotropic refinement
                  // along the face. If the situation
                  // described is not present, the data
                  // will be overwritten later on when we
                  // visit cells on finer levels, so no
                  // harm will be done.
                  if (face->has_children() &&
                      (cell->is_active() ||
                       GeometryInfo<dim>::face_refinement_case(
                         cell->refinement_case(), f) ==
                         RefinementCase<dim - 1>::isotropic_refinement))
                    {
                      for (unsigned int c = 0; c < face->n_children(); ++c)
                        adjacent_cells[2 * face->child(c)->index() + offset] =
                          cell;
                      if (face->child(0)->has_children())
                        {
                          adjacent_cells[2 * face->child(0)->child(0)->index() +
                                         offset] = cell;
                          adjacent_cells[2 * face->child(0)->child(1)->index() +
                                         offset] = cell;
                        }
                      if (face->child(1)->has_children())
                        {
                          adjacent_cells[2 * face->child(1)->child(0)->index() +
                                         offset] = cell;
                          adjacent_cells[2 * face->child(1)->child(1)->index() +
                                         offset] = cell;
                        }
                    } // if cell active and face refined
                }     // else -> dim==3
            }         // for all faces of all cells
 
        // now loop again over all cells and set the
        // corresponding neighbor cell. Note, that we
        // have to use the opposite of the
        // left_right_offset in this case as we want
        // the offset of the neighbor, not our own.
        for (const auto &cell : triangulation.cell_iterators())
          for (auto f : cell->face_indices())
            {
              const unsigned int offset =
                (cell->direction_flag() ?
                   left_right_offset[dim - 2][f][cell->face_orientation(f)] :
                   1 -
                     left_right_offset[dim - 2][f][cell->face_orientation(f)]);
              cell->set_neighbor(
                f, adjacent_cells[2 * cell->face(f)->index() + 1 - offset]);
            }
      }
 
 
      template <int dim, int spacedim>
      static void
      create_triangulation(const std::vector<Point<spacedim>> &vertices,
                           const std::vector<CellData<dim>>   &cells,
                           const SubCellData                  &subcelldata,
                           Triangulation<dim, spacedim>       &tria)
      {
        AssertThrow(vertices.size() > 0, ExcMessage("No vertices given"));
        AssertThrow(cells.size() > 0, ExcMessage("No cells given"));
 
        // Check that all cells have positive volume.
#ifndef _MSC_VER
        // TODO: The following code does not compile with MSVC. Find a way
        // around it
        if (dim == spacedim)
          for (unsigned int cell_no = 0; cell_no < cells.size(); ++cell_no)
            {
              // If we should check for distorted cells, then we permit them
              // to exist. If a cell has negative measure, then it must be
              // distorted (the converse is not necessarily true); hence
              // throw an exception if no such cells should exist.
              if (tria.check_for_distorted_cells)
                {
                  const double cell_measure = GridTools::cell_measure<spacedim>(
                    vertices,
                    ArrayView<const unsigned int>(cells[cell_no].vertices));
                  AssertThrow(cell_measure > 0, ExcGridHasInvalidCell(cell_no));
                }
            }
#endif
 
        // clear old content
        tria.levels.clear();
        tria.levels.push_back(
          std::make_unique<
            ::internal::TriangulationImplementation::TriaLevel>(dim));
 
        if (dim > 1)
          tria.faces = std::make_unique<
            ::internal::TriangulationImplementation::TriaFaces>(dim);
 
        // copy vertices
        tria.vertices = vertices;
        tria.vertices_used.assign(vertices.size(), true);
 
        // compute connectivity
        const auto connectivity   = build_connectivity<unsigned int>(cells);
        const unsigned int n_cell = cells.size();
 
        // TriaObjects: lines
        if (dim >= 2)
          {
            auto &lines_0 = tria.faces->lines; // data structure to be filled
 
            // get connectivity between quads and lines
            const auto        &crs     = connectivity.entity_to_entities(1, 0);
            const unsigned int n_lines = crs.ptr.size() - 1;
 
            // allocate memory
            reserve_space_(lines_0, n_lines);
 
            // loop over lines
            for (unsigned int line = 0; line < n_lines; ++line)
              for (unsigned int i = crs.ptr[line], j = 0; i < crs.ptr[line + 1];
                   ++i, ++j)
                lines_0.cells[line * GeometryInfo<1>::faces_per_cell + j] =
                  crs.col[i]; // set vertex indices
          }
 
        // TriaObjects: quads
        if (dim == 3)
          {
            auto &quads_0 = tria.faces->quads; // data structures to be filled
            auto &faces   = *tria.faces;
 
            // get connectivity between quads and lines
            const auto        &crs     = connectivity.entity_to_entities(2, 1);
            const unsigned int n_quads = crs.ptr.size() - 1;
 
            // allocate memory
            reserve_space_(quads_0, n_quads);
            reserve_space_(faces, 2 /*structdim*/, n_quads);
 
            // loop over all quads -> entity type, line indices/orientations
            for (unsigned int q = 0, k = 0; q < n_quads; ++q)
              {
                // set entity type of quads
                faces.set_quad_type(q, connectivity.entity_types(2)[q]);
 
                // loop over all its lines
                for (unsigned int i = crs.ptr[q], j = 0; i < crs.ptr[q + 1];
                     ++i, ++j, ++k)
                  {
                    // set line index
                    quads_0.cells[q * GeometryInfo<3>::lines_per_face + j] =
                      crs.col[i];
 
                    // set line orientations
                    const unsigned char combined_orientation =
                      connectivity.entity_orientations(1)
                        .get_combined_orientation(k);
                    // it doesn't make sense to set any flags except
                    // orientation for a line
                    Assert(
                      combined_orientation ==
                          ReferenceCell::default_combined_face_orientation() ||
                        combined_orientation ==
                          ReferenceCell::reversed_combined_line_orientation(),
                      ExcInternalError());
                    faces.quads_line_orientations
                      [q * GeometryInfo<3>::lines_per_face + j] =
                      combined_orientation ==
                      ReferenceCell::default_combined_face_orientation();
                  }
              }
          }
 
        // TriaObjects/TriaLevel: cell
        {
          auto &cells_0 = tria.levels[0]->cells; // data structure to be filled
          auto &level   = *tria.levels[0];
 
          // get connectivity between cells/faces and cells/cells
          const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
          const auto &nei = connectivity.entity_to_entities(dim, dim);
 
          // in 2d optional: since in in pure QUAD meshes same line
          // orientations can be guaranteed
          bool orientation_needed = false;
          if (dim == 3)
            orientation_needed = true;
          else if (dim == 2)
            {
              const auto &orientations = connectivity.entity_orientations(1);
              for (unsigned int i = 0; i < orientations.n_objects(); ++i)
                if (orientations.get_combined_orientation(i) !=
                    ReferenceCell::default_combined_face_orientation())
                  {
                    orientation_needed = true;
                    break;
                  }
            }
 
          // allocate memory
          reserve_space_(cells_0, n_cell);
          reserve_space_(level, spacedim, n_cell, orientation_needed);
 
          // loop over all cells
          for (unsigned int cell = 0; cell < n_cell; ++cell)
            {
              // set material ids
              cells_0.boundary_or_material_id[cell].material_id =
                cells[cell].material_id;
 
              // set manifold ids
              cells_0.manifold_id[cell] = cells[cell].manifold_id;
 
              // set entity types
              level.reference_cell[cell] = connectivity.entity_types(dim)[cell];
 
              // loop over faces
              for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
                   ++i, ++j)
                {
                  // set neighbor if not at boundary
                  if (nei.col[i] != static_cast<unsigned int>(-1))
                    level.neighbors[cell * GeometryInfo<dim>::faces_per_cell +
                                    j] = {0, nei.col[i]};
 
                  // set face indices
                  cells_0.cells[cell * GeometryInfo<dim>::faces_per_cell + j] =
                    crs.col[i];
 
                  // set face orientation if needed
                  if (orientation_needed)
                    {
                      level.face_orientations.set_combined_orientation(
                        cell * GeometryInfo<dim>::faces_per_cell + j,
                        connectivity.entity_orientations(dim - 1)
                          .get_combined_orientation(i));
                    }
                }
            }
        }
 
        // TriaFaces: boundary id of boundary faces
        if (dim > 1)
          {
            auto &bids_face = dim == 3 ?
                                tria.faces->quads.boundary_or_material_id :
                                tria.faces->lines.boundary_or_material_id;
 
            // count number of cells a face is belonging to
            std::vector<unsigned int> count(bids_face.size(), 0);
 
            // get connectivity between cells/faces
            const auto &crs = connectivity.entity_to_entities(dim, dim - 1);
 
            // count how many cells are adjacent to the same face
            for (unsigned int cell = 0; cell < cells.size(); ++cell)
              for (unsigned int i = crs.ptr[cell]; i < crs.ptr[cell + 1]; ++i)
                count[crs.col[i]]++;
 
            // loop over all faces
            for (unsigned int face = 0; face < count.size(); ++face)
              {
                if (count[face] != 1) // inner face
                  continue;
 
                // boundary faces ...
                bids_face[face].boundary_id = 0;
 
                if (dim != 3)
                  continue;
 
                // ... and the lines of quads in 3d
                const auto &crs = connectivity.entity_to_entities(2, 1);
                for (unsigned int i = crs.ptr[face]; i < crs.ptr[face + 1]; ++i)
                  tria.faces->lines.boundary_or_material_id[crs.col[i]]
                    .boundary_id = 0;
              }
          }
        else // 1d
          {
            static const unsigned int t_tba   = static_cast<unsigned int>(-1);
            static const unsigned int t_inner = static_cast<unsigned int>(-2);
 
            std::vector<unsigned int> type(vertices.size(), t_tba);
 
            const auto &crs = connectivity.entity_to_entities(1, 0);
 
            for (unsigned int cell = 0; cell < cells.size(); ++cell)
              for (unsigned int i = crs.ptr[cell], j = 0; i < crs.ptr[cell + 1];
                   ++i, ++j)
                if (type[crs.col[i]] != t_inner)
                  type[crs.col[i]] = type[crs.col[i]] == t_tba ? j : t_inner;
 
            for (unsigned int face = 0; face < type.size(); ++face)
              {
                // note: we also treat manifolds here!?
                (*tria.vertex_to_manifold_id_map_1d)[face] =
                  numbers::flat_manifold_id;
                if (type[face] != t_inner && type[face] != t_tba)
                  (*tria.vertex_to_boundary_id_map_1d)[face] = type[face];
              }
          }
 
        // SubCellData: line
        if (dim >= 2)
          process_subcelldata(connectivity.entity_to_entities(1, 0),
                              tria.faces->lines,
                              subcelldata.boundary_lines,
                              vertices);
 
        // SubCellData: quad
        if (dim == 3)
          process_subcelldata(connectivity.entity_to_entities(2, 0),
                              tria.faces->quads,
                              subcelldata.boundary_quads,
                              vertices);
      }
 
 
      template <int structdim, int spacedim, typename T>
      static void
      process_subcelldata(
        const CRS<T>                           &crs,
        TriaObjects                            &obj,
        const std::vector<CellData<structdim>> &boundary_objects_in,
        const std::vector<Point<spacedim>>     &vertex_locations)
      {
        AssertDimension(obj.structdim, structdim);
 
        if (boundary_objects_in.empty())
          return; // empty subcelldata -> nothing to do
 
        // pre-sort subcelldata
        auto boundary_objects = boundary_objects_in;
 
        // ... sort vertices
        for (auto &boundary_object : boundary_objects)
          std::sort(boundary_object.vertices.begin(),
                    boundary_object.vertices.end());
 
        // ... sort cells
        std::sort(boundary_objects.begin(),
                  boundary_objects.end(),
                  [](const auto &a, const auto &b) {
                    return a.vertices < b.vertices;
                  });
 
        unsigned int counter = 0;
 
        std::vector<unsigned int> key;
        key.reserve(GeometryInfo<structdim>::vertices_per_cell);
 
        for (unsigned int o = 0; o < obj.n_objects(); ++o)
          {
            auto &boundary_id = obj.boundary_or_material_id[o].boundary_id;
            auto &manifold_id = obj.manifold_id[o];
 
            // assert that object has not been visited yet and its value
            // has not been modified yet
            AssertThrow(boundary_id == 0 ||
                          boundary_id == numbers::internal_face_boundary_id,
                        ExcNotImplemented());
            AssertThrow(manifold_id == numbers::flat_manifold_id,
                        ExcNotImplemented());
 
            // create key
            key.assign(crs.col.data() + crs.ptr[o],
                       crs.col.data() + crs.ptr[o + 1]);
            std::sort(key.begin(), key.end());
 
            // is subcelldata provided? -> binary search
            const auto subcell_object =
              std::lower_bound(boundary_objects.begin(),
                               boundary_objects.end(),
                               key,
                               [&](const auto &cell, const auto &key) {
                                 return cell.vertices < key;
                               });
 
            // no subcelldata provided for this object
            if (subcell_object == boundary_objects.end() ||
                subcell_object->vertices != key)
              continue;
 
            ++counter;
 
            // set manifold id
            manifold_id = subcell_object->manifold_id;
 
            // set boundary id
            if (subcell_object->boundary_id !=
                numbers::internal_face_boundary_id)
              {
                (void)vertex_locations;
                AssertThrow(
                  boundary_id != numbers::internal_face_boundary_id,
                  ExcMessage(
                    "The input arguments for creating a triangulation "
                    "specified a boundary id for an internal face. This "
                    "is not allowed."
                    "\n\n"
                    "The object in question has vertex indices " +
                    [subcell_object]() {
                      std::string s;
                      for (const auto v : subcell_object->vertices)
                        s += std::to_string(v) + ',';
                      return s;
                    }() +
                    " which are located at positions " +
                    [vertex_locations, subcell_object]() {
                      std::ostringstream s;
                      for (const auto v : subcell_object->vertices)
                        s << '(' << vertex_locations[v] << ')';
                      return s.str();
                    }() +
                    "."));
                boundary_id = subcell_object->boundary_id;
              }
          }
 
        // make sure that all subcelldata entries have been processed
        // TODO: this is not guaranteed, why?
        // AssertDimension(counter, boundary_objects_in.size());
        (void)counter;
      }
 
 
 
      static void
      reserve_space_(TriaFaces         &faces,
                     const unsigned     structdim,
                     const unsigned int size)
      {
        const unsigned int dim = faces.dim;
 
        const unsigned int max_lines_per_face = 2 * structdim;
 
        if (dim == 3 && structdim == 2)
          {
            // quad entity types
            faces.quad_is_quadrilateral.assign(size, true);
 
            // quad line orientations
            faces.quads_line_orientations.assign(size * max_lines_per_face,
                                                 true);
          }
      }
 
 
 
      static void
      reserve_space_(TriaLevel         &level,
                     const unsigned int spacedim,
                     const unsigned int size,
                     const bool         orientation_needed)
      {
        const unsigned int dim = level.dim;
 
        const unsigned int max_faces_per_cell = 2 * dim;
 
        level.active_cell_indices.assign(size, numbers::invalid_unsigned_int);
        level.subdomain_ids.assign(size, 0);
        level.level_subdomain_ids.assign(size, 0);
 
        level.refine_flags.assign(size, 0u);
        level.refine_choice.assign(size, 0u);
        level.coarsen_flags.assign(size, false);
 
        level.parents.assign((size + 1) / 2, -1);
 
        if (dim == spacedim - 1)
          level.direction_flags.assign(size, true);
 
        level.neighbors.assign(size * max_faces_per_cell, {-1, -1});
 
        level.reference_cell.assign(size, ReferenceCells::Invalid);
 
        if (orientation_needed)
          level.face_orientations.reinit(size * max_faces_per_cell);
 
 
        level.global_active_cell_indices.assign(size,
                                                numbers::invalid_dof_index);
        level.global_level_cell_indices.assign(size,
                                               numbers::invalid_dof_index);
      }
 
 
 
      static void
      reserve_space_(TriaObjects &obj, const unsigned int size)
      {
        const unsigned int structdim = obj.structdim;
 
        const unsigned int max_children_per_cell = 1 << structdim;
        const unsigned int max_faces_per_cell    = 2 * structdim;
 
        obj.used.assign(size, true);
        obj.boundary_or_material_id.assign(
          size,
          internal::TriangulationImplementation::TriaObjects::
            BoundaryOrMaterialId());
        obj.manifold_id.assign(size, -1);
        obj.user_flags.assign(size, false);
        obj.user_data.resize(size);
 
        if (structdim > 1) // TODO: why?
          obj.refinement_cases.assign(size, 0);
 
        obj.children.assign(max_children_per_cell / 2 * size, -1);
 
        obj.cells.assign(max_faces_per_cell * size, -1);
 
        if (structdim <= 2)
          {
            obj.next_free_single               = size - 1;
            obj.next_free_pair                 = 0;
            obj.reverse_order_next_free_single = true;
          }
        else
          {
            obj.next_free_single = obj.next_free_pair = 0;
          }
      }
 
 
      template <int spacedim>
      static void
      delete_children(Triangulation<1, spacedim> &triangulation,
                      typename Triangulation<1, spacedim>::cell_iterator &cell,
                      std::vector<unsigned int> &,
                      std::vector<unsigned int> &)
      {
        const unsigned int dim = 1;
 
        // first we need to reset the
        // neighbor pointers of the
        // neighbors of this cell's
        // children to this cell. This is
        // different for one dimension,
        // since there neighbors can have a
        // refinement level differing from
        // that of this cell's children by
        // more than one level.
 
        Assert(!cell->child(0)->has_children() &&
                 !cell->child(1)->has_children(),
               ExcInternalError());
 
        // first do it for the cells to the
        // left
        if (cell->neighbor(0).state() == IteratorState::valid)
          if (cell->neighbor(0)->has_children())
            {
              typename Triangulation<dim, spacedim>::cell_iterator neighbor =
                cell->neighbor(0);
              Assert(neighbor->level() == cell->level(), ExcInternalError());
 
              // right child
              neighbor = neighbor->child(1);
              while (true)
                {
                  Assert(neighbor->neighbor(1) == cell->child(0),
                         ExcInternalError());
                  neighbor->set_neighbor(1, cell);
 
                  // move on to further
                  // children on the
                  // boundary between this
                  // cell and its neighbor
                  if (neighbor->has_children())
                    neighbor = neighbor->child(1);
                  else
                    break;
                }
            }
 
        // now do it for the cells to the
        // left
        if (cell->neighbor(1).state() == IteratorState::valid)
          if (cell->neighbor(1)->has_children())
            {
              typename Triangulation<dim, spacedim>::cell_iterator neighbor =
                cell->neighbor(1);
              Assert(neighbor->level() == cell->level(), ExcInternalError());
 
              // left child
              neighbor = neighbor->child(0);
              while (true)
                {
                  Assert(neighbor->neighbor(0) == cell->child(1),
                         ExcInternalError());
                  neighbor->set_neighbor(0, cell);
 
                  // move on to further
                  // children on the
                  // boundary between this
                  // cell and its neighbor
                  if (neighbor->has_children())
                    neighbor = neighbor->child(0);
                  else
                    break;
                }
            }
 
 
        // delete the vertex which will not
        // be needed anymore. This vertex
        // is the second of the first child
        triangulation.vertices_used[cell->child(0)->vertex_index(1)] = false;
 
        // invalidate children.  clear user
        // pointers, to avoid that they may
        // appear at unwanted places later
        // on...
        for (unsigned int child = 0; child < cell->n_children(); ++child)
          {
            cell->child(child)->clear_user_data();
            cell->child(child)->clear_user_flag();
            cell->child(child)->clear_used_flag();
          }
 
 
        // delete pointer to children
        cell->clear_children();
        cell->clear_user_flag();
      }
 
 
 
      template <int spacedim>
      static void
      delete_children(Triangulation<2, spacedim> &triangulation,
                      typename Triangulation<2, spacedim>::cell_iterator &cell,
                      std::vector<unsigned int> &line_cell_count,
                      std::vector<unsigned int> &)
      {
        const unsigned int        dim      = 2;
        const RefinementCase<dim> ref_case = cell->refinement_case();
 
        Assert(line_cell_count.size() == triangulation.n_raw_lines(),
               ExcInternalError());
 
        // vectors to hold all lines which
        // may be deleted
        std::vector<typename Triangulation<dim, spacedim>::line_iterator>
          lines_to_delete(0);
 
        lines_to_delete.reserve(4 * 2 + 4);
 
        // now we decrease the counters for
        // lines contained in the child
        // cells
        for (unsigned int c = 0; c < cell->n_children(); ++c)
          {
            typename Triangulation<dim, spacedim>::cell_iterator child =
              cell->child(c);
            for (unsigned int l = 0; l < GeometryInfo<dim>::lines_per_cell; ++l)
              --line_cell_count[child->line_index(l)];
          }
 
 
        // delete the vertex which will not
        // be needed anymore. This vertex
        // is the second of the second line
        // of the first child, if the cell
        // is refined with cut_xy, else there
        // is no inner vertex.
        // additionally delete unneeded inner
        // lines
        if (ref_case == RefinementCase<dim>::cut_xy)
          {
            triangulation
              .vertices_used[cell->child(0)->line(1)->vertex_index(1)] = false;
 
            lines_to_delete.push_back(cell->child(0)->line(1));
            lines_to_delete.push_back(cell->child(0)->line(3));
            lines_to_delete.push_back(cell->child(3)->line(0));
            lines_to_delete.push_back(cell->child(3)->line(2));
          }
        else
          {
            unsigned int inner_face_no =
              ref_case == RefinementCase<dim>::cut_x ? 1 : 3;
 
            // the inner line will not be
            // used any more
            lines_to_delete.push_back(cell->child(0)->line(inner_face_no));
          }
 
        // invalidate children
        for (unsigned int child = 0; child < cell->n_children(); ++child)
          {
            cell->child(child)->clear_user_data();
            cell->child(child)->clear_user_flag();
            cell->child(child)->clear_used_flag();
          }
 
 
        // delete pointer to children
        cell->clear_children();
        cell->clear_refinement_case();
        cell->clear_user_flag();
 
        // look at the refinement of outer
        // lines. if nobody needs those
        // anymore we can add them to the
        // list of lines to be deleted.
        for (unsigned int line_no = 0;
             line_no < GeometryInfo<dim>::lines_per_cell;
             ++line_no)
          {
            typename Triangulation<dim, spacedim>::line_iterator line =
              cell->line(line_no);
 
            if (line->has_children())
              {
                // if one of the cell counters is
                // zero, the other has to be as well
 
                Assert((line_cell_count[line->child_index(0)] == 0 &&
                        line_cell_count[line->child_index(1)] == 0) ||
                         (line_cell_count[line->child_index(0)] > 0 &&
                          line_cell_count[line->child_index(1)] > 0),
                       ExcInternalError());
 
                if (line_cell_count[line->child_index(0)] == 0)
                  {
                    for (unsigned int c = 0; c < 2; ++c)
                      Assert(!line->child(c)->has_children(),
                             ExcInternalError());
 
                    // we may delete the line's
                    // children and the middle vertex
                    // as no cell references them
                    // anymore
                    triangulation
                      .vertices_used[line->child(0)->vertex_index(1)] = false;
 
                    lines_to_delete.push_back(line->child(0));
                    lines_to_delete.push_back(line->child(1));
 
                    line->clear_children();
                  }
              }
          }
 
        // finally, delete unneeded lines
 
        // clear user pointers, to avoid that
        // they may appear at unwanted places
        // later on...
        // same for user flags, then finally
        // delete the lines
        typename std::vector<
          typename Triangulation<dim, spacedim>::line_iterator>::iterator
          line    = lines_to_delete.begin(),
          endline = lines_to_delete.end();
        for (; line != endline; ++line)
          {
            (*line)->clear_user_data();
            (*line)->clear_user_flag();
            (*line)->clear_used_flag();
          }
      }
 
 
 
      template <int spacedim>
      static void
      delete_children(Triangulation<3, spacedim> &triangulation,
                      typename Triangulation<3, spacedim>::cell_iterator &cell,
                      std::vector<unsigned int> &line_cell_count,
                      std::vector<unsigned int> &quad_cell_count)
      {
        const unsigned int dim = 3;
 
        Assert(line_cell_count.size() == triangulation.n_raw_lines(),
               ExcInternalError());
        Assert(quad_cell_count.size() == triangulation.n_raw_quads(),
               ExcInternalError());
 
        // first of all, we store the RefineCase of
        // this cell
        const RefinementCase<dim> ref_case = cell->refinement_case();
        // vectors to hold all lines and quads which
        // may be deleted
        std::vector<typename Triangulation<dim, spacedim>::line_iterator>
          lines_to_delete(0);
        std::vector<typename Triangulation<dim, spacedim>::quad_iterator>
          quads_to_delete(0);
 
        lines_to_delete.reserve(12 * 2 + 6 * 4 + 6);
        quads_to_delete.reserve(6 * 4 + 12);
 
        // now we decrease the counters for lines and
        // quads contained in the child cells
        for (unsigned int c = 0; c < cell->n_children(); ++c)
          {
            typename Triangulation<dim, spacedim>::cell_iterator child =
              cell->child(c);
            const auto line_indices = TriaAccessorImplementation::
              Implementation::get_line_indices_of_cell(*child);
            for (const unsigned int l : cell->line_indices())
              --line_cell_count[line_indices[l]];
            for (auto f : GeometryInfo<dim>::face_indices())
              --quad_cell_count[child->quad_index(f)];
          }
 
        //-------------------------------------
        // delete interior quads and lines and the
        // interior vertex, depending on the
        // refinement case of the cell
        //
        // for append quads and lines: only append
        // them to the list of objects to be deleted
 
        switch (ref_case)
          {
            case RefinementCase<dim>::cut_x:
              quads_to_delete.push_back(cell->child(0)->face(1));
              break;
            case RefinementCase<dim>::cut_y:
              quads_to_delete.push_back(cell->child(0)->face(3));
              break;
            case RefinementCase<dim>::cut_z:
              quads_to_delete.push_back(cell->child(0)->face(5));
              break;
            case RefinementCase<dim>::cut_xy:
              quads_to_delete.push_back(cell->child(0)->face(1));
              quads_to_delete.push_back(cell->child(0)->face(3));
              quads_to_delete.push_back(cell->child(3)->face(0));
              quads_to_delete.push_back(cell->child(3)->face(2));
 
              lines_to_delete.push_back(cell->child(0)->line(11));
              break;
            case RefinementCase<dim>::cut_xz:
              quads_to_delete.push_back(cell->child(0)->face(1));
              quads_to_delete.push_back(cell->child(0)->face(5));
              quads_to_delete.push_back(cell->child(3)->face(0));
              quads_to_delete.push_back(cell->child(3)->face(4));
 
              lines_to_delete.push_back(cell->child(0)->line(5));
              break;
            case RefinementCase<dim>::cut_yz:
              quads_to_delete.push_back(cell->child(0)->face(3));
              quads_to_delete.push_back(cell->child(0)->face(5));
              quads_to_delete.push_back(cell->child(3)->face(2));
              quads_to_delete.push_back(cell->child(3)->face(4));
 
              lines_to_delete.push_back(cell->child(0)->line(7));
              break;
            case RefinementCase<dim>::cut_xyz:
              quads_to_delete.push_back(cell->child(0)->face(1));
              quads_to_delete.push_back(cell->child(2)->face(1));
              quads_to_delete.push_back(cell->child(4)->face(1));
              quads_to_delete.push_back(cell->child(6)->face(1));
 
              quads_to_delete.push_back(cell->child(0)->face(3));
              quads_to_delete.push_back(cell->child(1)->face(3));
              quads_to_delete.push_back(cell->child(4)->face(3));
              quads_to_delete.push_back(cell->child(5)->face(3));
 
              quads_to_delete.push_back(cell->child(0)->face(5));
              quads_to_delete.push_back(cell->child(1)->face(5));
              quads_to_delete.push_back(cell->child(2)->face(5));
              quads_to_delete.push_back(cell->child(3)->face(5));
 
              lines_to_delete.push_back(cell->child(0)->line(5));
              lines_to_delete.push_back(cell->child(0)->line(7));
              lines_to_delete.push_back(cell->child(0)->line(11));
              lines_to_delete.push_back(cell->child(7)->line(0));
              lines_to_delete.push_back(cell->child(7)->line(2));
              lines_to_delete.push_back(cell->child(7)->line(8));
              // delete the vertex which will not
              // be needed anymore. This vertex
              // is the vertex at the heart of
              // this cell, which is the sixth of
              // the first child
              triangulation.vertices_used[cell->child(0)->vertex_index(7)] =
                false;
              break;
            default:
              // only remaining case is
              // no_refinement, thus an error
              DEAL_II_ASSERT_UNREACHABLE();
              break;
          }
 
 
        // invalidate children
        for (unsigned int child = 0; child < cell->n_children(); ++child)
          {
            cell->child(child)->clear_user_data();
            cell->child(child)->clear_user_flag();
 
            for (auto f : GeometryInfo<dim>::face_indices())
              // set flags denoting deviations from standard orientation of
              // faces back to initialization values
              cell->child(child)->set_combined_face_orientation(
                f, ReferenceCell::default_combined_face_orientation());
 
            cell->child(child)->clear_used_flag();
          }
 
 
        // delete pointer to children
        cell->clear_children();
        cell->clear_refinement_case();
        cell->clear_user_flag();
 
        // so far we only looked at inner quads,
        // lines and vertices. Now we have to
        // consider outer ones as well. here, we have
        // to check, whether there are other cells
        // still needing these objects. otherwise we
        // can delete them. first for quads (and
        // their inner lines).
 
        for (const unsigned int quad_no : GeometryInfo<dim>::face_indices())
          {
            typename Triangulation<dim, spacedim>::quad_iterator quad =
              cell->face(quad_no);
 
            Assert(
              (GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) &&
               quad->has_children()) ||
                GeometryInfo<dim>::face_refinement_case(ref_case, quad_no) ==
                  RefinementCase<dim - 1>::no_refinement,
              ExcInternalError());
 
            switch (quad->refinement_case())
              {
                case RefinementCase<dim - 1>::no_refinement:
                  // nothing to do as the quad
                  // is not refined
                  break;
                case RefinementCase<dim - 1>::cut_x:
                case RefinementCase<dim - 1>::cut_y:
                  {
                    // if one of the cell counters is
                    // zero, the other has to be as
                    // well
                    Assert((quad_cell_count[quad->child_index(0)] == 0 &&
                            quad_cell_count[quad->child_index(1)] == 0) ||
                             (quad_cell_count[quad->child_index(0)] > 0 &&
                              quad_cell_count[quad->child_index(1)] > 0),
                           ExcInternalError());
                    // it might be, that the quad is
                    // refined twice anisotropically,
                    // first check, whether we may
                    // delete possible grand_children
                    unsigned int deleted_grandchildren       = 0;
                    unsigned int number_of_child_refinements = 0;
 
                    for (unsigned int c = 0; c < 2; ++c)
                      if (quad->child(c)->has_children())
                        {
                          ++number_of_child_refinements;
                          // if one of the cell counters is
                          // zero, the other has to be as
                          // well
                          Assert(
                            (quad_cell_count[quad->child(c)->child_index(0)] ==
                               0 &&
                             quad_cell_count[quad->child(c)->child_index(1)] ==
                               0) ||
                              (quad_cell_count[quad->child(c)->child_index(0)] >
                                 0 &&
                               quad_cell_count[quad->child(c)->child_index(1)] >
                                 0),
                            ExcInternalError());
                          if (quad_cell_count[quad->child(c)->child_index(0)] ==
                              0)
                            {
                              // Assert, that the two
                              // anisotropic
                              // refinements add up to
                              // isotropic refinement
                              Assert(quad->refinement_case() +
                                         quad->child(c)->refinement_case() ==
                                       RefinementCase<dim>::cut_xy,
                                     ExcInternalError());
                              // we may delete the
                              // quad's children and
                              // the inner line as no
                              // cell references them
                              // anymore
                              quads_to_delete.push_back(
                                quad->child(c)->child(0));
                              quads_to_delete.push_back(
                                quad->child(c)->child(1));
                              if (quad->child(c)->refinement_case() ==
                                  RefinementCase<2>::cut_x)
                                lines_to_delete.push_back(
                                  quad->child(c)->child(0)->line(1));
                              else
                                lines_to_delete.push_back(
                                  quad->child(c)->child(0)->line(3));
                              quad->child(c)->clear_children();
                              quad->child(c)->clear_refinement_case();
                              ++deleted_grandchildren;
                            }
                        }
                    // if no grandchildren are left, we
                    // may as well delete the
                    // refinement of the inner line
                    // between our children and the
                    // corresponding vertex
                    if (number_of_child_refinements > 0 &&
                        deleted_grandchildren == number_of_child_refinements)
                      {
                        typename Triangulation<dim, spacedim>::line_iterator
                          middle_line;
                        if (quad->refinement_case() == RefinementCase<2>::cut_x)
                          middle_line = quad->child(0)->line(1);
                        else
                          middle_line = quad->child(0)->line(3);
 
                        lines_to_delete.push_back(middle_line->child(0));
                        lines_to_delete.push_back(middle_line->child(1));
                        triangulation
                          .vertices_used[middle_vertex_index<dim, spacedim>(
                            middle_line)] = false;
                        middle_line->clear_children();
                      }
 
                    // now consider the direct children
                    // of the given quad
                    if (quad_cell_count[quad->child_index(0)] == 0)
                      {
                        // we may delete the quad's
                        // children and the inner line
                        // as no cell references them
                        // anymore
                        quads_to_delete.push_back(quad->child(0));
                        quads_to_delete.push_back(quad->child(1));
                        if (quad->refinement_case() == RefinementCase<2>::cut_x)
                          lines_to_delete.push_back(quad->child(0)->line(1));
                        else
                          lines_to_delete.push_back(quad->child(0)->line(3));
 
                        // if the counters just dropped
                        // to zero, otherwise the
                        // children would have been
                        // deleted earlier, then this
                        // cell's children must have
                        // contained the anisotropic
                        // quad children. thus, if
                        // those have again anisotropic
                        // children, which are in
                        // effect isotropic children of
                        // the original quad, those are
                        // still needed by a
                        // neighboring cell and we
                        // cannot delete them. instead,
                        // we have to reset this quad's
                        // refine case to isotropic and
                        // set the children
                        // accordingly.
                        if (quad->child(0)->has_children())
                          if (quad->refinement_case() ==
                              RefinementCase<2>::cut_x)
                            {
                              // now evereything is
                              // quite complicated. we
                              // have the children
                              // numbered according to
                              //
                              // *---*---*
                              // |n+1|m+1|
                              // *---*---*
                              // | n | m |
                              // *---*---*
                              //
                              // from the original
                              // anisotropic
                              // refinement. we have to
                              // reorder them as
                              //
                              // *---*---*
                              // | m |m+1|
                              // *---*---*
                              // | n |n+1|
                              // *---*---*
                              //
                              // for isotropic refinement.
                              //
                              // this is a bit ugly, of
                              // course: loop over all
                              // cells on all levels
                              // and look for faces n+1
                              // (switch_1) and m
                              // (switch_2).
                              const typename Triangulation<dim, spacedim>::
                                quad_iterator switch_1 =
                                                quad->child(0)->child(1),
                                              switch_2 =
                                                quad->child(1)->child(0);
 
                              Assert(!switch_1->has_children(),
                                     ExcInternalError());
                              Assert(!switch_2->has_children(),
                                     ExcInternalError());
 
                              const int switch_1_index = switch_1->index();
                              const int switch_2_index = switch_2->index();
                              for (unsigned int l = 0;
                                   l < triangulation.levels.size();
                                   ++l)
                                for (unsigned int h = 0;
                                     h <
                                     triangulation.levels[l]->cells.n_objects();
                                     ++h)
                                  for (const unsigned int q :
                                       GeometryInfo<dim>::face_indices())
                                    {
                                      const int index =
                                        triangulation.levels[l]
                                          ->cells.get_bounding_object_indices(
                                            h)[q];
                                      if (index == switch_1_index)
                                        triangulation.levels[l]
                                          ->cells.get_bounding_object_indices(
                                            h)[q] = switch_2_index;
                                      else if (index == switch_2_index)
                                        triangulation.levels[l]
                                          ->cells.get_bounding_object_indices(
                                            h)[q] = switch_1_index;
                                    }
                              // now we have to copy
                              // all information of the
                              // two quads
                              const int switch_1_lines[4] = {
                                static_cast<signed int>(
                                  switch_1->line_index(0)),
                                static_cast<signed int>(
                                  switch_1->line_index(1)),
                                static_cast<signed int>(
                                  switch_1->line_index(2)),
                                static_cast<signed int>(
                                  switch_1->line_index(3))};
                              const bool switch_1_line_orientations[4] = {
                                switch_1->line_orientation(0),
                                switch_1->line_orientation(1),
                                switch_1->line_orientation(2),
                                switch_1->line_orientation(3)};
                              const types::boundary_id switch_1_boundary_id =
                                switch_1->boundary_id();
                              const unsigned int switch_1_user_index =
                                switch_1->user_index();
                              const bool switch_1_user_flag =
                                switch_1->user_flag_set();
 
                              switch_1->set_bounding_object_indices(
                                {switch_2->line_index(0),
                                 switch_2->line_index(1),
                                 switch_2->line_index(2),
                                 switch_2->line_index(3)});
                              switch_1->set_line_orientation(
                                0, switch_2->line_orientation(0));
                              switch_1->set_line_orientation(
                                1, switch_2->line_orientation(1));
                              switch_1->set_line_orientation(
                                2, switch_2->line_orientation(2));
                              switch_1->set_line_orientation(
                                3, switch_2->line_orientation(3));
                              switch_1->set_boundary_id_internal(
                                switch_2->boundary_id());
                              switch_1->set_manifold_id(
                                switch_2->manifold_id());
                              switch_1->set_user_index(switch_2->user_index());
                              if (switch_2->user_flag_set())
                                switch_1->set_user_flag();
                              else
                                switch_1->clear_user_flag();
 
                              switch_2->set_bounding_object_indices(
                                {switch_1_lines[0],
                                 switch_1_lines[1],
                                 switch_1_lines[2],
                                 switch_1_lines[3]});
                              switch_2->set_line_orientation(
                                0, switch_1_line_orientations[0]);
                              switch_2->set_line_orientation(
                                1, switch_1_line_orientations[1]);
                              switch_2->set_line_orientation(
                                2, switch_1_line_orientations[2]);
                              switch_2->set_line_orientation(
                                3, switch_1_line_orientations[3]);
                              switch_2->set_boundary_id_internal(
                                switch_1_boundary_id);
                              switch_2->set_manifold_id(
                                switch_1->manifold_id());
                              switch_2->set_user_index(switch_1_user_index);
                              if (switch_1_user_flag)
                                switch_2->set_user_flag();
                              else
                                switch_2->clear_user_flag();
 
                              const unsigned int child_0 =
                                quad->child(0)->child_index(0);
                              const unsigned int child_2 =
                                quad->child(1)->child_index(0);
                              quad->clear_children();
                              quad->clear_refinement_case();
                              quad->set_refinement_case(
                                RefinementCase<2>::cut_xy);
                              quad->set_children(0, child_0);
                              quad->set_children(2, child_2);
                              std::swap(quad_cell_count[child_0 + 1],
                                        quad_cell_count[child_2]);
                            }
                          else
                            {
                              // the face was refined
                              // with cut_y, thus the
                              // children are already
                              // in correct order. we
                              // only have to set them
                              // correctly, deleting
                              // the indirection of two
                              // anisotropic refinement
                              // and going directly
                              // from the quad to
                              // isotropic children
                              const unsigned int child_0 =
                                quad->child(0)->child_index(0);
                              const unsigned int child_2 =
                                quad->child(1)->child_index(0);
                              quad->clear_children();
                              quad->clear_refinement_case();
                              quad->set_refinement_case(
                                RefinementCase<2>::cut_xy);
                              quad->set_children(0, child_0);
                              quad->set_children(2, child_2);
                            }
                        else
                          {
                            quad->clear_children();
                            quad->clear_refinement_case();
                          }
                      }
                    break;
                  }
                case RefinementCase<dim - 1>::cut_xy:
                  {
                    // if one of the cell counters is
                    // zero, the others have to be as
                    // well
 
                    Assert((quad_cell_count[quad->child_index(0)] == 0 &&
                            quad_cell_count[quad->child_index(1)] == 0 &&
                            quad_cell_count[quad->child_index(2)] == 0 &&
                            quad_cell_count[quad->child_index(3)] == 0) ||
                             (quad_cell_count[quad->child_index(0)] > 0 &&
                              quad_cell_count[quad->child_index(1)] > 0 &&
                              quad_cell_count[quad->child_index(2)] > 0 &&
                              quad_cell_count[quad->child_index(3)] > 0),
                           ExcInternalError());
 
                    if (quad_cell_count[quad->child_index(0)] == 0)
                      {
                        // we may delete the quad's
                        // children, the inner lines
                        // and the middle vertex as no
                        // cell references them anymore
                        lines_to_delete.push_back(quad->child(0)->line(1));
                        lines_to_delete.push_back(quad->child(3)->line(0));
                        lines_to_delete.push_back(quad->child(0)->line(3));
                        lines_to_delete.push_back(quad->child(3)->line(2));
 
                        for (unsigned int child = 0; child < quad->n_children();
                             ++child)
                          quads_to_delete.push_back(quad->child(child));
 
                        triangulation
                          .vertices_used[quad->child(0)->vertex_index(3)] =
                          false;
 
                        quad->clear_children();
                        quad->clear_refinement_case();
                      }
                  }
                  break;
 
                default:
                  DEAL_II_ASSERT_UNREACHABLE();
                  break;
              }
          }
 
        // now we repeat a similar procedure
        // for the outer lines of this cell.
 
        // if in debug mode: check that each
        // of the lines for which we consider
        // deleting the children in fact has
        // children (the bits/coarsening_3d
        // test tripped over this initially)
        for (unsigned int line_no = 0;
             line_no < GeometryInfo<dim>::lines_per_cell;
             ++line_no)
          {
            typename Triangulation<dim, spacedim>::line_iterator line =
              cell->line(line_no);
 
            Assert(
              (GeometryInfo<dim>::line_refinement_case(ref_case, line_no) &&
               line->has_children()) ||
                GeometryInfo<dim>::line_refinement_case(ref_case, line_no) ==
                  RefinementCase<1>::no_refinement,
              ExcInternalError());
 
            if (line->has_children())
              {
                // if one of the cell counters is
                // zero, the other has to be as well
 
                Assert((line_cell_count[line->child_index(0)] == 0 &&
                        line_cell_count[line->child_index(1)] == 0) ||
                         (line_cell_count[line->child_index(0)] > 0 &&
                          line_cell_count[line->child_index(1)] > 0),
                       ExcInternalError());
 
                if (line_cell_count[line->child_index(0)] == 0)
                  {
                    for (unsigned int c = 0; c < 2; ++c)
                      Assert(!line->child(c)->has_children(),
                             ExcInternalError());
 
                    // we may delete the line's
                    // children and the middle vertex
                    // as no cell references them
                    // anymore
                    triangulation
                      .vertices_used[line->child(0)->vertex_index(1)] = false;
 
                    lines_to_delete.push_back(line->child(0));
                    lines_to_delete.push_back(line->child(1));
 
                    line->clear_children();
                  }
              }
          }
 
        // finally, delete unneeded quads and lines
 
        // clear user pointers, to avoid that
        // they may appear at unwanted places
        // later on...
        // same for user flags, then finally
        // delete the quads and lines
        typename std::vector<
          typename Triangulation<dim, spacedim>::line_iterator>::iterator
          line    = lines_to_delete.begin(),
          endline = lines_to_delete.end();
        for (; line != endline; ++line)
          {
            (*line)->clear_user_data();
            (*line)->clear_user_flag();
            (*line)->clear_used_flag();
          }
 
        typename std::vector<
          typename Triangulation<dim, spacedim>::quad_iterator>::iterator
          quad    = quads_to_delete.begin(),
          endquad = quads_to_delete.end();
        for (; quad != endquad; ++quad)
          {
            (*quad)->clear_user_data();
            (*quad)->clear_children();
            (*quad)->clear_refinement_case();
            (*quad)->clear_user_flag();
            (*quad)->clear_used_flag();
          }
      }
 
 
      template <int spacedim>
      static void
      create_children(
        Triangulation<2, spacedim> &triangulation,
        unsigned int               &next_unused_vertex,
        typename Triangulation<2, spacedim>::raw_line_iterator
          &next_unused_line,
        typename Triangulation<2, spacedim>::raw_cell_iterator
          &next_unused_cell,
        const typename Triangulation<2, spacedim>::cell_iterator &cell)
      {
        const unsigned int dim = 2;
        // clear refinement flag
        const RefinementCase<dim> ref_case = cell->refine_flag_set();
        cell->clear_refine_flag();
 
        /* For the refinement process: since we go the levels up from the
           lowest, there are (unlike above) only two possibilities: a neighbor
           cell is on the same level or one level up (in both cases, it may or
           may not be refined later on, but we don't care here).
 
           First:
           Set up an array of the 3x3 vertices, which are distributed on the
           cell (the array consists of indices into the @p{vertices} std::vector
 
           2--7--3
           |  |  |
           4--8--5
           |  |  |
           0--6--1
 
           note: in case of cut_x or cut_y not all these vertices are needed for
           the new cells
 
           Second:
           Set up an array of the new lines (the array consists of iterator
           pointers into the lines arrays)
 
           .-6-.-7-.         The directions are:  .->-.->-.
           1   9   3                              ^   ^   ^
           .-10.11-.                             .->-.->-.
           0   8   2                              ^   ^   ^
           .-4-.-5-.                              .->-.->-.
 
           cut_x:
           .-4-.-5-.
           |   |   |
           0   6   1
           |   |   |
           .-2-.-3-.
 
           cut_y:
           .---5---.
           1       3
           .---6---.
           0       2
           .---4---.
 
 
           Third:
           Set up an array of neighbors:
 
           6  7
           .--.--.
           1|  |  |3
           .--.--.
           0|  |  |2
           .--.--.
           4   5
 
           We need this array for two reasons: first to get the lines which will
           bound the four subcells (if the neighboring cell is refined, these
           lines already exist), and second to update neighborship information.
           Since if a neighbor is not refined, its neighborship record only
           points to the present, unrefined, cell rather than the children we
           are presently creating, we only need the neighborship information
           if the neighbor cells are refined. In all other cases, we store
           the unrefined neighbor address
 
           We also need for every neighbor (if refined) which number among its
           neighbors the present (unrefined) cell has, since that number is to
           be replaced and because that also is the number of the subline which
           will be the interface between that neighbor and the to be created
           cell. We will store this number (between 0 and 3) in the field
           @p{neighbors_neighbor}.
 
           It would be sufficient to use the children of the common line to the
           neighbor, if we only wanted to get the new sublines and the new
           vertex, but because we need to update the neighborship information of
           the two refined subcells of the neighbor, we need to search these
           anyway.
 
           Convention:
           The created children are numbered like this:
 
           .--.--.
           |2 . 3|
           .--.--.
           |0 | 1|
           .--.--.
        */
        // collect the indices of the eight surrounding vertices
        //   2--7--3
        //   |  |  |
        //   4--8--5
        //   |  |  |
        //   0--6--1
        int new_vertices[9];
        for (unsigned int vertex_no = 0; vertex_no < 4; ++vertex_no)
          new_vertices[vertex_no] = cell->vertex_index(vertex_no);
        for (unsigned int line_no = 0; line_no < 4; ++line_no)
          if (cell->line(line_no)->has_children())
            new_vertices[4 + line_no] =
              cell->line(line_no)->child(0)->vertex_index(1);
 
        if (ref_case == RefinementCase<dim>::cut_xy)
          {
            // find the next
            // unused vertex and
            // allocate it for
            // the new vertex we
            // need here
            while (triangulation.vertices_used[next_unused_vertex] == true)
              ++next_unused_vertex;
            Assert(next_unused_vertex < triangulation.vertices.size(),
                   ExcMessage(
                     "Internal error: During refinement, the triangulation "
                     "wants to access an element of the 'vertices' array "
                     "but it turns out that the array is not large enough."));
            triangulation.vertices_used[next_unused_vertex] = true;
 
            new_vertices[8] = next_unused_vertex;
 
            // determine middle vertex by transfinite interpolation to be
            // consistent with what happens to quads in a
            // Triangulation<3,3> when they are refined
            triangulation.vertices[next_unused_vertex] =
              cell->center(true, true);
          }
 
 
        // Now the lines:
        typename Triangulation<dim, spacedim>::raw_line_iterator new_lines[12];
        unsigned int                                             lmin = 8;
        unsigned int                                             lmax = 12;
        if (ref_case != RefinementCase<dim>::cut_xy)
          {
            lmin = 6;
            lmax = 7;
          }
 
        for (unsigned int l = lmin; l < lmax; ++l)
          {
            while (next_unused_line->used() == true)
              ++next_unused_line;
            new_lines[l] = next_unused_line;
            ++next_unused_line;
 
            AssertIsNotUsed(new_lines[l]);
          }
 
        if (ref_case == RefinementCase<dim>::cut_xy)
          {
            //   .-6-.-7-.
            //   1   9   3
            //   .-10.11-.
            //   0   8   2
            //   .-4-.-5-.
 
            // lines 0-7 already exist, create only the four interior
            // lines 8-11
            unsigned int l = 0;
            for (const unsigned int face_no : GeometryInfo<dim>::face_indices())
              for (unsigned int c = 0; c < 2; ++c, ++l)
                new_lines[l] = cell->line(face_no)->child(c);
            Assert(l == 8, ExcInternalError());
 
            new_lines[8]->set_bounding_object_indices(
              {new_vertices[6], new_vertices[8]});
            new_lines[9]->set_bounding_object_indices(
              {new_vertices[8], new_vertices[7]});
            new_lines[10]->set_bounding_object_indices(
              {new_vertices[4], new_vertices[8]});
            new_lines[11]->set_bounding_object_indices(
              {new_vertices[8], new_vertices[5]});
          }
        else if (ref_case == RefinementCase<dim>::cut_x)
          {
            //   .-4-.-5-.
            //   |   |   |
            //   0   6   1
            //   |   |   |
            //   .-2-.-3-.
            new_lines[0] = cell->line(0);
            new_lines[1] = cell->line(1);
            new_lines[2] = cell->line(2)->child(0);
            new_lines[3] = cell->line(2)->child(1);
            new_lines[4] = cell->line(3)->child(0);
            new_lines[5] = cell->line(3)->child(1);
            new_lines[6]->set_bounding_object_indices(
              {new_vertices[6], new_vertices[7]});
          }
        else
          {
            Assert(ref_case == RefinementCase<dim>::cut_y, ExcInternalError());
            //   .---5---.
            //   1       3
            //   .---6---.
            //   0       2
            //   .---4---.
            new_lines[0] = cell->line(0)->child(0);
            new_lines[1] = cell->line(0)->child(1);
            new_lines[2] = cell->line(1)->child(0);
            new_lines[3] = cell->line(1)->child(1);
            new_lines[4] = cell->line(2);
            new_lines[5] = cell->line(3);
            new_lines[6]->set_bounding_object_indices(
              {new_vertices[4], new_vertices[5]});
          }
 
        for (unsigned int l = lmin; l < lmax; ++l)
          {
            new_lines[l]->set_used_flag();
            new_lines[l]->clear_user_flag();
            new_lines[l]->clear_user_data();
            new_lines[l]->clear_children();
            // interior line
            new_lines[l]->set_boundary_id_internal(
              numbers::internal_face_boundary_id);
            new_lines[l]->set_manifold_id(cell->manifold_id());
          }
 
        // Now add the four (two)
        // new cells!
        typename Triangulation<dim, spacedim>::raw_cell_iterator
          subcells[GeometryInfo<dim>::max_children_per_cell];
        while (next_unused_cell->used() == true)
          ++next_unused_cell;
 
        const unsigned int n_children = GeometryInfo<dim>::n_children(ref_case);
        for (unsigned int i = 0; i < n_children; ++i)
          {
            AssertIsNotUsed(next_unused_cell);
            subcells[i] = next_unused_cell;
            ++next_unused_cell;
            if (i % 2 == 1 && i < n_children - 1)
              while (next_unused_cell->used() == true)
                ++next_unused_cell;
          }
 
        if (ref_case == RefinementCase<dim>::cut_xy)
          {
            // children:
            //   .--.--.
            //   |2 . 3|
            //   .--.--.
            //   |0 | 1|
            //   .--.--.
            // lines:
            //   .-6-.-7-.
            //   1   9   3
            //   .-10.11-.
            //   0   8   2
            //   .-4-.-5-.
            subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
                                                      new_lines[8]->index(),
                                                      new_lines[4]->index(),
                                                      new_lines[10]->index()});
            subcells[1]->set_bounding_object_indices({new_lines[8]->index(),
                                                      new_lines[2]->index(),
                                                      new_lines[5]->index(),
                                                      new_lines[11]->index()});
            subcells[2]->set_bounding_object_indices({new_lines[1]->index(),
                                                      new_lines[9]->index(),
                                                      new_lines[10]->index(),
                                                      new_lines[6]->index()});
            subcells[3]->set_bounding_object_indices({new_lines[9]->index(),
                                                      new_lines[3]->index(),
                                                      new_lines[11]->index(),
                                                      new_lines[7]->index()});
          }
        else if (ref_case == RefinementCase<dim>::cut_x)
          {
            // children:
            //   .--.--.
            //   |  .  |
            //   .0 . 1.
            //   |  |  |
            //   .--.--.
            // lines:
            //   .-4-.-5-.
            //   |   |   |
            //   0   6   1
            //   |   |   |
            //   .-2-.-3-.
            subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
                                                      new_lines[6]->index(),
                                                      new_lines[2]->index(),
                                                      new_lines[4]->index()});
            subcells[1]->set_bounding_object_indices({new_lines[6]->index(),
                                                      new_lines[1]->index(),
                                                      new_lines[3]->index(),
                                                      new_lines[5]->index()});
          }
        else
          {
            Assert(ref_case == RefinementCase<dim>::cut_y, ExcInternalError());
            // children:
            //   .-----.
            //   |  1  |
            //   .-----.
            //   |  0  |
            //   .-----.
            // lines:
            //   .---5---.
            //   1       3
            //   .---6---.
            //   0       2
            //   .---4---.
            subcells[0]->set_bounding_object_indices({new_lines[0]->index(),
                                                      new_lines[2]->index(),
                                                      new_lines[4]->index(),
                                                      new_lines[6]->index()});
            subcells[1]->set_bounding_object_indices({new_lines[1]->index(),
                                                      new_lines[3]->index(),
                                                      new_lines[6]->index(),
                                                      new_lines[5]->index()});
          }
 
        types::subdomain_id subdomainid = cell->subdomain_id();
 
        for (unsigned int i = 0; i < n_children; ++i)
          {
            subcells[i]->set_used_flag();
            subcells[i]->clear_refine_flag();
            subcells[i]->clear_user_flag();
            subcells[i]->clear_user_data();
            subcells[i]->clear_children();
            // inherit material properties
            subcells[i]->set_material_id(cell->material_id());
            subcells[i]->set_manifold_id(cell->manifold_id());
            subcells[i]->set_subdomain_id(subdomainid);
 
            if (i % 2 == 0)
              subcells[i]->set_parent(cell->index());
          }
 
 
 
        // set child index for even children i=0,2 (0)
        for (unsigned int i = 0; i < n_children / 2; ++i)
          cell->set_children(2 * i, subcells[2 * i]->index());
        // set the refine case
        cell->set_refinement_case(ref_case);
 
        // note that the
        // refinement flag was
        // already cleared at the
        // beginning of this function
 
        if (dim == spacedim - 1)
          for (unsigned int c = 0; c < n_children; ++c)
            cell->child(c)->set_direction_flag(cell->direction_flag());
      }
 
 
 
      template <int dim, int spacedim>
      static typename Triangulation<dim, spacedim>::DistortedCellList
      execute_refinement_isotropic(Triangulation<dim, spacedim> &triangulation,
                                   const bool check_for_distorted_cells)
      {
        AssertDimension(dim, 2);
 
        // Check whether a new level is needed. We have to check for
        // this on the highest level only
        for (const auto &cell : triangulation.active_cell_iterators_on_level(
               triangulation.levels.size() - 1))
          if (cell->refine_flag_set())
            {
              triangulation.levels.push_back(
                std::make_unique<
                  internal::TriangulationImplementation::TriaLevel>(dim));
              break;
            }
 
        for (typename Triangulation<dim, spacedim>::line_iterator line =
               triangulation.begin_line();
             line != triangulation.end_line();
             ++line)
          {
            line->clear_user_flag();
            line->clear_user_data();
          }
 
        unsigned int n_single_lines   = 0;
        unsigned int n_lines_in_pairs = 0;
        unsigned int needed_vertices  = 0;
 
        for (int level = triangulation.levels.size() - 2; level >= 0; --level)
          {
            // count number of flagged cells on this level and compute
            // how many new vertices and new lines will be needed
            unsigned int needed_cells = 0;
 
            for (const auto &cell :
                 triangulation.active_cell_iterators_on_level(level))
              if (cell->refine_flag_set())
                {
                  if (cell->reference_cell() == ReferenceCells::Triangle)
                    {
                      needed_cells += 4;
                      needed_vertices += 0;
                      n_single_lines += 3;
                    }
                  else if (cell->reference_cell() ==
                           ReferenceCells::Quadrilateral)
                    {
                      needed_cells += 4;
                      needed_vertices += 1;
                      n_single_lines += 4;
                    }
                  else
                    {
                      AssertThrow(false, ExcNotImplemented());
                    }
 
                  for (const auto line_no : cell->face_indices())
                    {
                      auto line = cell->line(line_no);
                      if (line->has_children() == false)
                        line->set_user_flag();
                    }
                }
 
 
            const unsigned int used_cells =
              std::count(triangulation.levels[level + 1]->cells.used.begin(),
                         triangulation.levels[level + 1]->cells.used.end(),
                         true);
 
 
            reserve_space(*triangulation.levels[level + 1],
                          used_cells + needed_cells,
                          2,
                          spacedim);
 
            reserve_space(triangulation.levels[level + 1]->cells,
                          needed_cells,
                          0);
          }
 
        for (auto line = triangulation.begin_line();
             line != triangulation.end_line();
             ++line)
          if (line->user_flag_set())
            {
              Assert(line->has_children() == false, ExcInternalError());
              n_lines_in_pairs += 2;
              needed_vertices += 1;
            }
 
        reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
 
        needed_vertices += std::count(triangulation.vertices_used.begin(),
                                      triangulation.vertices_used.end(),
                                      true);
 
        if (needed_vertices > triangulation.vertices.size())
          {
            triangulation.vertices.resize(needed_vertices, Point<spacedim>());
            triangulation.vertices_used.resize(needed_vertices, false);
          }
 
        unsigned int next_unused_vertex = 0;
 
        {
          typename Triangulation<dim, spacedim>::active_line_iterator
            line = triangulation.begin_active_line(),
            endl = triangulation.end_line();
          typename Triangulation<dim, spacedim>::raw_line_iterator
            next_unused_line = triangulation.begin_raw_line();
 
          for (; line != endl; ++line)
            if (line->user_flag_set())
              {
                // This line needs to be refined. Find the next unused vertex
                // and set it appropriately
                while (triangulation.vertices_used[next_unused_vertex] == true)
                  ++next_unused_vertex;
                Assert(next_unused_vertex < triangulation.vertices.size(),
                       ExcMessage(
                         "Internal error: During refinement, the triangulation "
                         "wants to access an element of the 'vertices' array "
                         "but it turns out that the array is not large "
                         "enough."));
                triangulation.vertices_used[next_unused_vertex] = true;
 
                triangulation.vertices[next_unused_vertex] = line->center(true);
 
                bool pair_found = false;
                (void)pair_found;
                for (; next_unused_line != endl; ++next_unused_line)
                  if (!next_unused_line->used() &&
                      !(++next_unused_line)->used())
                    {
                      --next_unused_line;
                      pair_found = true;
                      break;
                    }
                Assert(pair_found, ExcInternalError());
 
                line->set_children(0, next_unused_line->index());
 
                const typename Triangulation<dim, spacedim>::raw_line_iterator
                  children[2] = {next_unused_line, ++next_unused_line};
 
                AssertIsNotUsed(children[0]);
                AssertIsNotUsed(children[1]);
 
                children[0]->set_bounding_object_indices(
                  {line->vertex_index(0), next_unused_vertex});
                children[1]->set_bounding_object_indices(
                  {next_unused_vertex, line->vertex_index(1)});
 
                for (auto &child : children)
                  {
                    child->set_used_flag();
                    child->clear_children();
                    child->clear_user_data();
                    child->clear_user_flag();
                    child->set_boundary_id_internal(line->boundary_id());
                    child->set_manifold_id(line->manifold_id());
                    // Line orientation is relative to the cell it is on so
                    // those cannot be set at this point.
                  }
 
                line->clear_user_flag();
              }
        }
 
        reserve_space(triangulation.faces->lines, 0, n_single_lines);
 
        typename Triangulation<dim, spacedim>::DistortedCellList
          cells_with_distorted_children;
 
        typename Triangulation<dim, spacedim>::raw_line_iterator
          next_unused_line = triangulation.begin_raw_line();
 
        const auto create_children = [](auto         &triangulation,
                                        unsigned int &next_unused_vertex,
                                        auto         &next_unused_line,
                                        auto         &next_unused_cell,
                                        const auto   &cell) {
          const auto ref_case = cell->refine_flag_set();
          cell->clear_refine_flag();
 
          unsigned int n_new_vertices = 0;
 
          if (cell->reference_cell() == ReferenceCells::Triangle)
            n_new_vertices = 6;
          else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
            n_new_vertices = 9;
          else
            AssertThrow(false, ExcNotImplemented());
 
          std::vector<unsigned int> new_vertices(n_new_vertices,
                                                 numbers::invalid_unsigned_int);
          for (unsigned int vertex_no = 0; vertex_no < cell->n_vertices();
               ++vertex_no)
            new_vertices[vertex_no] = cell->vertex_index(vertex_no);
          for (unsigned int line_no = 0; line_no < cell->n_lines(); ++line_no)
            if (cell->line(line_no)->has_children())
              new_vertices[cell->n_vertices() + line_no] =
                cell->line(line_no)->child(0)->vertex_index(1);
 
          if (cell->reference_cell() == ReferenceCells::Quadrilateral)
            {
              while (triangulation.vertices_used[next_unused_vertex] == true)
                ++next_unused_vertex;
              Assert(
                next_unused_vertex < triangulation.vertices.size(),
                ExcMessage(
                  "Internal error: During refinement, the triangulation wants "
                  "to access an element of the 'vertices' array but it turns "
                  "out that the array is not large enough."));
              triangulation.vertices_used[next_unused_vertex] = true;
 
              new_vertices[8] = next_unused_vertex;
 
              triangulation.vertices[next_unused_vertex] =
                cell->center(true, true);
            }
 
          std::array<typename Triangulation<dim, spacedim>::raw_line_iterator,
                     12>
                                        new_lines;
          std::array<unsigned char, 12> inherited_orientations;
          inherited_orientations.fill(
            ReferenceCell::default_combined_face_orientation());
          unsigned int lmin = 0;
          unsigned int lmax = 0;
 
          if (cell->reference_cell() == ReferenceCells::Triangle)
            {
              lmin = 6;
              lmax = 9;
              // For triangles, the innermost faces are always reversed for the
              // first three children and are in the standard orientation for
              // the last one.
              std::fill(inherited_orientations.begin() + lmin,
                        inherited_orientations.begin() + lmax,
                        ReferenceCell::reversed_combined_line_orientation());
            }
          else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
            {
              lmin = 8;
              lmax = 12;
            }
          else
            {
              AssertThrow(false, ExcNotImplemented());
            }
 
          for (unsigned int l = lmin; l < lmax; ++l)
            {
              while (next_unused_line->used() == true)
                ++next_unused_line;
              new_lines[l] = next_unused_line;
              ++next_unused_line;
 
              AssertIsNotUsed(new_lines[l]);
            }
 
          // set up lines which have parents:
          for (const unsigned int face_no : cell->face_indices())
            {
              // Check the face (line) orientation to ensure that the (six or
              // eight) outer lines in new_lines are indexed in the default
              // orientation. This way we can index into this array in the
              // without special casing orientations (e.g., quadrilateral child
              // 3 will always have lines 9, 3, 11, 7) when setting child lines.
              const unsigned char combined_orientation =
                cell->combined_face_orientation(face_no);
              Assert(combined_orientation ==
                         ReferenceCell::default_combined_face_orientation() ||
                       combined_orientation ==
                         ReferenceCell::reversed_combined_line_orientation(),
                     ExcInternalError());
              for (unsigned int c = 0; c < 2; ++c)
                {
                  new_lines[2 * face_no + c] = cell->line(face_no)->child(c);
                  inherited_orientations[2 * face_no + c] =
                    cell->combined_face_orientation(face_no);
                }
              if (combined_orientation ==
                  ReferenceCell::reversed_combined_line_orientation())
                std::swap(new_lines[2 * face_no], new_lines[2 * face_no + 1]);
            }
 
          // set up lines which do not have parents:
          if (cell->reference_cell() == ReferenceCells::Triangle)
            {
              new_lines[6]->set_bounding_object_indices(
                {new_vertices[3], new_vertices[4]});
              new_lines[7]->set_bounding_object_indices(
                {new_vertices[4], new_vertices[5]});
              new_lines[8]->set_bounding_object_indices(
                {new_vertices[5], new_vertices[3]});
            }
          else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
            {
              new_lines[8]->set_bounding_object_indices(
                {new_vertices[6], new_vertices[8]});
              new_lines[9]->set_bounding_object_indices(
                {new_vertices[8], new_vertices[7]});
              new_lines[10]->set_bounding_object_indices(
                {new_vertices[4], new_vertices[8]});
              new_lines[11]->set_bounding_object_indices(
                {new_vertices[8], new_vertices[5]});
            }
          else
            {
              AssertThrow(false, ExcNotImplemented());
            }
 
          for (unsigned int l = lmin; l < lmax; ++l)
            {
              new_lines[l]->set_used_flag();
              new_lines[l]->clear_user_flag();
              new_lines[l]->clear_user_data();
              new_lines[l]->clear_children();
              // new lines are always internal.
              new_lines[l]->set_boundary_id_internal(
                numbers::internal_face_boundary_id);
              new_lines[l]->set_manifold_id(cell->manifold_id());
            }
 
          typename Triangulation<dim, spacedim>::raw_cell_iterator
            subcells[GeometryInfo<dim>::max_children_per_cell];
          while (next_unused_cell->used() == true)
            ++next_unused_cell;
 
          unsigned int n_children = 0;
          if (cell->reference_cell() == ReferenceCells::Triangle)
            n_children = 4;
          else if (cell->reference_cell() == ReferenceCells::Quadrilateral)
            n_children = 4;
          else
            AssertThrow(false, ExcNotImplemented());
 
          for (unsigned int i = 0; i < n_children; ++i)
            {
              AssertIsNotUsed(next_unused_cell);
              subcells[i] = next_unused_cell;
              ++next_unused_cell;
              if (i % 2 == 1 && i < n_children - 1)
                while (next_unused_cell->used() == true)
                  ++next_unused_cell;
            }
 
          // Assign lines to child cells:
          constexpr unsigned int X = numbers::invalid_unsigned_int;
          static constexpr ::ndarray<unsigned int, 4, 4> tri_child_lines =
            {{{{0, 8, 5, X}}, {{1, 2, 6, X}}, {{7, 3, 4, X}}, {{6, 7, 8, X}}}};
          static constexpr ::ndarray<unsigned int, 4, 4>
            quad_child_lines = {{{{0, 8, 4, 10}},
                                 {{8, 2, 5, 11}},
                                 {{1, 9, 10, 6}},
                                 {{9, 3, 11, 7}}}};
          // Here and below we assume that child cells have the same reference
          // cell type as the parent.
          const auto &child_lines =
            cell->reference_cell() == ReferenceCells::Triangle ?
              tri_child_lines :
              quad_child_lines;
          for (unsigned int i = 0; i < n_children; ++i)
            {
              if (cell->reference_cell() == ReferenceCells::Triangle)
                subcells[i]->set_bounding_object_indices(
                  {new_lines[child_lines[i][0]]->index(),
                   new_lines[child_lines[i][1]]->index(),
                   new_lines[child_lines[i][2]]->index()});
              else
                subcells[i]->set_bounding_object_indices(
                  {new_lines[child_lines[i][0]]->index(),
                   new_lines[child_lines[i][1]]->index(),
                   new_lines[child_lines[i][2]]->index(),
                   new_lines[child_lines[i][3]]->index()});
 
              subcells[i]->set_used_flag();
              subcells[i]->clear_refine_flag();
              subcells[i]->clear_user_flag();
              subcells[i]->clear_user_data();
              subcells[i]->clear_children();
              // inherit material properties
              subcells[i]->set_material_id(cell->material_id());
              subcells[i]->set_manifold_id(cell->manifold_id());
              subcells[i]->set_subdomain_id(cell->subdomain_id());
 
              triangulation.levels[subcells[i]->level()]
                ->reference_cell[subcells[i]->index()] = cell->reference_cell();
 
              // Finally, now that children are marked as used, we can set
              // orientation flags:
              for (unsigned int face_no : cell->face_indices())
                subcells[i]->set_combined_face_orientation(
                  face_no, inherited_orientations[child_lines[i][face_no]]);
 
              if (i % 2 == 0)
                subcells[i]->set_parent(cell->index());
            }
 
          // Unlike the same lines on other children, the innermost triangle's
          // faces are all in the default orientation:
          if (cell->reference_cell() == ReferenceCells::Triangle)
            for (unsigned int face_no : cell->face_indices())
              subcells[3]->set_combined_face_orientation(
                face_no, ReferenceCell::default_combined_face_orientation());
 
          for (unsigned int i = 0; i < n_children / 2; ++i)
            cell->set_children(2 * i, subcells[2 * i]->index());
 
          cell->set_refinement_case(ref_case);
 
          if (dim == spacedim - 1)
            for (unsigned int c = 0; c < n_children; ++c)
              cell->child(c)->set_direction_flag(cell->direction_flag());
        };
 
        for (int level = 0;
             level < static_cast<int>(triangulation.levels.size()) - 1;
             ++level)
          {
            typename Triangulation<dim, spacedim>::raw_cell_iterator
              next_unused_cell = triangulation.begin_raw(level + 1);
 
            for (const auto &cell :
                 triangulation.active_cell_iterators_on_level(level))
              if (cell->refine_flag_set())
                {
                  create_children(triangulation,
                                  next_unused_vertex,
                                  next_unused_line,
                                  next_unused_cell,
                                  cell);
 
                  if (cell->reference_cell() == ReferenceCells::Quadrilateral &&
                      check_for_distorted_cells &&
                      has_distorted_children<dim, spacedim>(cell))
                    cells_with_distorted_children.distorted_cells.push_back(
                      cell);
 
                  triangulation.signals.post_refinement_on_cell(cell);
                }
          }
 
        return cells_with_distorted_children;
      }
 
 
 
      template <int spacedim>
      static typename Triangulation<1, spacedim>::DistortedCellList
      execute_refinement(Triangulation<1, spacedim> &triangulation,
                         const bool /*check_for_distorted_cells*/)
      {
        const unsigned int dim = 1;
 
        // Check whether a new level is needed. We have to check for
        // this on the highest level only
        for (const auto &cell : triangulation.active_cell_iterators_on_level(
               triangulation.levels.size() - 1))
          if (cell->refine_flag_set())
            {
              triangulation.levels.push_back(
                std::make_unique<
                  internal::TriangulationImplementation::TriaLevel>(dim));
              break;
            }
 
 
        // check how much space is needed on every level. We need not
        // check the highest level since either - on the highest level
        // no cells are flagged for refinement - there are, but
        // prepare_refinement added another empty level
        unsigned int needed_vertices = 0;
        for (int level = triangulation.levels.size() - 2; level >= 0; --level)
          {
            // count number of flagged
            // cells on this level
            unsigned int flagged_cells = 0;
 
            for (const auto &acell :
                 triangulation.active_cell_iterators_on_level(level))
              if (acell->refine_flag_set())
                ++flagged_cells;
 
            // count number of used cells
            // on the next higher level
            const unsigned int used_cells =
              std::count(triangulation.levels[level + 1]->cells.used.begin(),
                         triangulation.levels[level + 1]->cells.used.end(),
                         true);
 
            // reserve space for the used_cells cells already existing
            // on the next higher level as well as for the
            // 2*flagged_cells that will be created on that level
            reserve_space(*triangulation.levels[level + 1],
                          used_cells + GeometryInfo<1>::max_children_per_cell *
                                         flagged_cells,
                          1,
                          spacedim);
            // reserve space for 2*flagged_cells new lines on the next
            // higher level
            reserve_space(triangulation.levels[level + 1]->cells,
                          GeometryInfo<1>::max_children_per_cell *
                            flagged_cells,
                          0);
 
            needed_vertices += flagged_cells;
          }
 
        // add to needed vertices how many
        // vertices are already in use
        needed_vertices += std::count(triangulation.vertices_used.begin(),
                                      triangulation.vertices_used.end(),
                                      true);
        // if we need more vertices: create them, if not: leave the
        // array as is, since shrinking is not really possible because
        // some of the vertices at the end may be in use
        if (needed_vertices > triangulation.vertices.size())
          {
            triangulation.vertices.resize(needed_vertices, Point<spacedim>());
            triangulation.vertices_used.resize(needed_vertices, false);
          }
 
 
        // Do REFINEMENT on every level; exclude highest level as
        // above
 
        // index of next unused vertex
        unsigned int next_unused_vertex = 0;
 
        for (int level = triangulation.levels.size() - 2; level >= 0; --level)
          {
            typename Triangulation<dim, spacedim>::raw_cell_iterator
              next_unused_cell = triangulation.begin_raw(level + 1);
 
            for (const auto &cell :
                 triangulation.active_cell_iterators_on_level(level))
              if (cell->refine_flag_set())
                {
                  // clear refinement flag
                  cell->clear_refine_flag();
 
                  // search for next unused
                  // vertex
                  while (triangulation.vertices_used[next_unused_vertex] ==
                         true)
                    ++next_unused_vertex;
                  Assert(
                    next_unused_vertex < triangulation.vertices.size(),
                    ExcMessage(
                      "Internal error: During refinement, the triangulation "
                      "wants to access an element of the 'vertices' array "
                      "but it turns out that the array is not large enough."));
 
                  // Now we always ask the cell itself where to put
                  // the new point. The cell in turn will query the
                  // manifold object internally.
                  triangulation.vertices[next_unused_vertex] =
                    cell->center(true);
 
                  triangulation.vertices_used[next_unused_vertex] = true;
 
                  // search for next two unused cell (++ takes care of
                  // the end of the vector)
                  typename Triangulation<dim, spacedim>::raw_cell_iterator
                    first_child,
                    second_child;
                  while (next_unused_cell->used() == true)
                    ++next_unused_cell;
                  first_child = next_unused_cell;
                  first_child->set_used_flag();
                  first_child->clear_user_data();
                  ++next_unused_cell;
                  AssertIsNotUsed(next_unused_cell);
                  second_child = next_unused_cell;
                  second_child->set_used_flag();
                  second_child->clear_user_data();
 
                  types::subdomain_id subdomainid = cell->subdomain_id();
 
                  // insert first child
                  cell->set_children(0, first_child->index());
                  first_child->clear_children();
                  first_child->set_bounding_object_indices(
                    {cell->vertex_index(0), next_unused_vertex});
                  first_child->set_material_id(cell->material_id());
                  first_child->set_manifold_id(cell->manifold_id());
                  first_child->set_subdomain_id(subdomainid);
                  if (dim == spacedim - 1)
                    first_child->set_direction_flag(cell->direction_flag());
 
                  first_child->set_parent(cell->index());
 
                  // Set manifold id of the right face. Only do this
                  // on the first child.
                  first_child->face(1)->set_manifold_id(cell->manifold_id());
 
                  // reset neighborship info (refer to
                  // internal::TriangulationImplementation::TriaLevel<0> for
                  // details)
                  first_child->set_neighbor(1, second_child);
                  if (cell->neighbor(0).state() != IteratorState::valid)
                    first_child->set_neighbor(0, cell->neighbor(0));
                  else if (cell->neighbor(0)->is_active())
                    {
                      // since the neighbors level is always <=level,
                      // if the cell is active, then there are no
                      // cells to the left which may want to know
                      // about this new child cell.
                      Assert(cell->neighbor(0)->level() <= cell->level(),
                             ExcInternalError());
                      first_child->set_neighbor(0, cell->neighbor(0));
                    }
                  else
                    // left neighbor is refined
                    {
                      // set neighbor to cell on same level
                      const unsigned int nbnb = cell->neighbor_of_neighbor(0);
                      first_child->set_neighbor(0,
                                                cell->neighbor(0)->child(nbnb));
 
                      // reset neighbor info of all right descendant
                      // of the left neighbor of cell
                      typename Triangulation<dim, spacedim>::cell_iterator
                        left_neighbor = cell->neighbor(0);
                      while (left_neighbor->has_children())
                        {
                          left_neighbor = left_neighbor->child(nbnb);
                          left_neighbor->set_neighbor(nbnb, first_child);
                        }
                    }
 
                  // insert second child
                  second_child->clear_children();
                  second_child->set_bounding_object_indices(
                    {next_unused_vertex, cell->vertex_index(1)});
                  second_child->set_neighbor(0, first_child);
                  second_child->set_material_id(cell->material_id());
                  second_child->set_manifold_id(cell->manifold_id());
                  second_child->set_subdomain_id(subdomainid);
                  if (dim == spacedim - 1)
                    second_child->set_direction_flag(cell->direction_flag());
 
                  if (cell->neighbor(1).state() != IteratorState::valid)
                    second_child->set_neighbor(1, cell->neighbor(1));
                  else if (cell->neighbor(1)->is_active())
                    {
                      Assert(cell->neighbor(1)->level() <= cell->level(),
                             ExcInternalError());
                      second_child->set_neighbor(1, cell->neighbor(1));
                    }
                  else
                    // right neighbor is refined same as above
                    {
                      const unsigned int nbnb = cell->neighbor_of_neighbor(1);
                      second_child->set_neighbor(
                        1, cell->neighbor(1)->child(nbnb));
 
                      typename Triangulation<dim, spacedim>::cell_iterator
                        right_neighbor = cell->neighbor(1);
                      while (right_neighbor->has_children())
                        {
                          right_neighbor = right_neighbor->child(nbnb);
                          right_neighbor->set_neighbor(nbnb, second_child);
                        }
                    }
                  // inform all listeners that cell refinement is done
                  triangulation.signals.post_refinement_on_cell(cell);
                }
          }
 
        // in 1d, we can not have distorted children unless the parent
        // was already distorted (that is because we don't use
        // boundary information for 1d triangulations). so return an
        // empty list
        return typename Triangulation<1, spacedim>::DistortedCellList();
      }
 
 
      template <int spacedim>
      static typename Triangulation<2, spacedim>::DistortedCellList
      execute_refinement(Triangulation<2, spacedim> &triangulation,
                         const bool                  check_for_distorted_cells)
      {
        const unsigned int dim = 2;
 
        // First check whether we can get away with isotropic refinement, or
        // whether we need to run through the full anisotropic algorithm
        {
          bool do_isotropic_refinement = true;
          for (const auto &cell : triangulation.active_cell_iterators())
            if (cell->refine_flag_set() == RefinementCase<dim>::cut_x ||
                cell->refine_flag_set() == RefinementCase<dim>::cut_y)
              {
                do_isotropic_refinement = false;
                break;
              }
 
          if (do_isotropic_refinement)
            return execute_refinement_isotropic(triangulation,
                                                check_for_distorted_cells);
        }
 
        // If we get here, we are doing anisotropic refinement.
 
        // Check whether a new level is needed. We have to check for
        // this on the highest level only
        for (const auto &cell : triangulation.active_cell_iterators_on_level(
               triangulation.levels.size() - 1))
          if (cell->refine_flag_set())
            {
              triangulation.levels.push_back(
                std::make_unique<
                  internal::TriangulationImplementation::TriaLevel>(dim));
              break;
            }
 
        // TODO[WB]: we clear user flags and pointers of lines; we're going
        // to use them to flag which lines need refinement
        for (typename Triangulation<dim, spacedim>::line_iterator line =
               triangulation.begin_line();
             line != triangulation.end_line();
             ++line)
          {
            line->clear_user_flag();
            line->clear_user_data();
          }
        // running over all cells and lines count the number
        // n_single_lines of lines which can be stored as single
        // lines, e.g. inner lines
        unsigned int n_single_lines = 0;
 
        // New lines to be created: number lines which are stored in
        // pairs (the children of lines must be stored in pairs)
        unsigned int n_lines_in_pairs = 0;
 
        // check how much space is needed on every level. We need not
        // check the highest level since either - on the highest level
        // no cells are flagged for refinement - there are, but
        // prepare_refinement added another empty level
        unsigned int needed_vertices = 0;
        for (int level = triangulation.levels.size() - 2; level >= 0; --level)
          {
            // count number of flagged cells on this level and compute
            // how many new vertices and new lines will be needed
            unsigned int needed_cells = 0;
 
            for (const auto &cell :
                 triangulation.active_cell_iterators_on_level(level))
              if (cell->refine_flag_set())
                {
                  if (cell->refine_flag_set() == RefinementCase<dim>::cut_xy)
                    {
                      needed_cells += 4;
 
                      // new vertex at center of cell is needed in any
                      // case
                      ++needed_vertices;
 
                      // the four inner lines can be stored as singles
                      n_single_lines += 4;
                    }
                  else // cut_x || cut_y
                    {
                      // set the flag showing that anisotropic
                      // refinement is used for at least one cell
                      triangulation.anisotropic_refinement = true;
 
                      needed_cells += 2;
                      // no vertex at center
 
                      // the inner line can be stored as single
                      n_single_lines += 1;
                    }
 
                  // mark all faces (lines) for refinement; checking
                  // locally whether the neighbor would also like to
                  // refine them is rather difficult for lines so we
                  // only flag them and after visiting all cells, we
                  // decide which lines need refinement;
                  for (const unsigned int line_no :
                       GeometryInfo<dim>::face_indices())
                    {
                      if (GeometryInfo<dim>::face_refinement_case(
                            cell->refine_flag_set(), line_no) ==
                          RefinementCase<1>::cut_x)
                        {
                          typename Triangulation<dim, spacedim>::line_iterator
                            line = cell->line(line_no);
                          if (line->has_children() == false)
                            line->set_user_flag();
                        }
                    }
                }
 
 
            // count number of used cells on the next higher level
            const unsigned int used_cells =
              std::count(triangulation.levels[level + 1]->cells.used.begin(),
                         triangulation.levels[level + 1]->cells.used.end(),
                         true);
 
 
            // reserve space for the used_cells cells already existing
            // on the next higher level as well as for the
            // needed_cells that will be created on that level
            reserve_space(*triangulation.levels[level + 1],
                          used_cells + needed_cells,
                          2,
                          spacedim);
 
            // reserve space for needed_cells new quads on the next
            // higher level
            reserve_space(triangulation.levels[level + 1]->cells,
                          needed_cells,
                          0);
          }
 
        // now count the lines which were flagged for refinement
        for (typename Triangulation<dim, spacedim>::line_iterator line =
               triangulation.begin_line();
             line != triangulation.end_line();
             ++line)
          if (line->user_flag_set())
            {
              Assert(line->has_children() == false, ExcInternalError());
              n_lines_in_pairs += 2;
              needed_vertices += 1;
            }
        // reserve space for n_lines_in_pairs new lines.  note, that
        // we can't reserve space for the single lines here as well,
        // as all the space reserved for lines in pairs would be
        // counted as unused and we would end up with too little space
        // to store all lines. memory reservation for n_single_lines
        // can only be done AFTER we refined the lines of the current
        // cells
        reserve_space(triangulation.faces->lines, n_lines_in_pairs, 0);
 
        // add to needed vertices how many vertices are already in use
        needed_vertices += std::count(triangulation.vertices_used.begin(),
                                      triangulation.vertices_used.end(),
                                      true);
        // if we need more vertices: create them, if not: leave the
        // array as is, since shrinking is not really possible because
        // some of the vertices at the end may be in use
        if (needed_vertices > triangulation.vertices.size())
          {
            triangulation.vertices.resize(needed_vertices, Point<spacedim>());
            triangulation.vertices_used.resize(needed_vertices, false);
          }
 
 
        // Do REFINEMENT on every level; exclude highest level as
        // above
 
        //  index of next unused vertex
        unsigned int next_unused_vertex = 0;
 
        // first the refinement of lines.  children are stored
        // pairwise
        {
          // only active objects can be refined further
          typename Triangulation<dim, spacedim>::active_line_iterator
            line = triangulation.begin_active_line(),
            endl = triangulation.end_line();
          typename Triangulation<dim, spacedim>::raw_line_iterator
            next_unused_line = triangulation.begin_raw_line();
 
          for (; line != endl; ++line)
            if (line->user_flag_set())
              {
                // this line needs to be refined
 
                // find the next unused vertex and set it
                // appropriately
                while (triangulation.vertices_used[next_unused_vertex] == true)
                  ++next_unused_vertex;
                Assert(
                  next_unused_vertex < triangulation.vertices.size(),
                  ExcMessage(
                    "Internal error: During refinement, the triangulation wants to access an element of the 'vertices' array but it turns out that the array is not large enough."));
                triangulation.vertices_used[next_unused_vertex] = true;
 
                triangulation.vertices[next_unused_vertex] = line->center(true);
 
                // now that we created the right point, make up the
                // two child lines.  To this end, find a pair of
                // unused lines
                bool pair_found = false;
                (void)pair_found;
                for (; next_unused_line != endl; ++next_unused_line)
                  if (!next_unused_line->used() &&
                      !(++next_unused_line)->used())
                    {
                      // go back to the first of the two unused
                      // lines
                      --next_unused_line;
                      pair_found = true;
                      break;
                    }
                Assert(pair_found, ExcInternalError());
 
                // there are now two consecutive unused lines, such
                // that the children of a line will be consecutive.
                // then set the child pointer of the present line
                line->set_children(0, next_unused_line->index());
 
                // set the two new lines
                const typename Triangulation<dim, spacedim>::raw_line_iterator
                  children[2] = {next_unused_line, ++next_unused_line};
                // some tests; if any of the iterators should be
                // invalid, then already dereferencing will fail
                AssertIsNotUsed(children[0]);
                AssertIsNotUsed(children[1]);
 
                children[0]->set_bounding_object_indices(
                  {line->vertex_index(0), next_unused_vertex});
                children[1]->set_bounding_object_indices(
                  {next_unused_vertex, line->vertex_index(1)});
 
                children[0]->set_used_flag();
                children[1]->set_used_flag();
                children[0]->clear_children();
                children[1]->clear_children();
                children[0]->clear_user_data();
                children[1]->clear_user_data();
                children[0]->clear_user_flag();
                children[1]->clear_user_flag();
 
 
                children[0]->set_boundary_id_internal(line->boundary_id());
                children[1]->set_boundary_id_internal(line->boundary_id());
 
                children[0]->set_manifold_id(line->manifold_id());
                children[1]->set_manifold_id(line->manifold_id());
 
                // finally clear flag indicating the need for
                // refinement
                line->clear_user_flag();
              }
        }
 
 
        // Now set up the new cells
 
        // reserve space for inner lines (can be stored as single
        // lines)
        reserve_space(triangulation.faces->lines, 0, n_single_lines);
 
        typename Triangulation<2, spacedim>::DistortedCellList
          cells_with_distorted_children;
 
        // reset next_unused_line, as now also single empty places in
        // the vector can be used
        typename Triangulation<dim, spacedim>::raw_line_iterator
          next_unused_line = triangulation.begin_raw_line();
 
        for (int level = 0;
             level < static_cast<int>(triangulation.levels.size()) - 1;
             ++level)
          {
            typename Triangulation<dim, spacedim>::raw_cell_iterator
              next_unused_cell = triangulation.begin_raw(level + 1);
 
            for (const auto &cell :
                 triangulation.active_cell_iterators_on_level(level))
              if (cell->refine_flag_set())
                {
                  // actually set up the children and update neighbor
                  // information
                  create_children(triangulation,
                                  next_unused_vertex,
                                  next_unused_line,
                                  next_unused_cell,
                                  cell);
 
                  if (check_for_distorted_cells &&
                      has_distorted_children<dim, spacedim>(cell))
                    cells_with_distorted_children.distorted_cells.push_back(
                      cell);
                  // inform all listeners that cell refinement is done
                  triangulation.signals.post_refinement_on_cell(cell);
                }
          }
 
        return cells_with_distorted_children;
      }
 
 
      template <int spacedim>
      static typename Triangulation<3, spacedim>::DistortedCellList
      execute_refinement_isotropic(Triangulation<3, spacedim> &triangulation,
                                   const bool check_for_distorted_cells)
      {
        static const int          dim = 3;
        static const unsigned int X   = numbers::invalid_unsigned_int;
        using raw_line_iterator =
          typename Triangulation<dim, spacedim>::raw_line_iterator;
        using raw_quad_iterator =
          typename Triangulation<dim, spacedim>::raw_quad_iterator;
 
        Assert(spacedim == 3, ExcNotImplemented());
 
        Assert(triangulation.vertices.size() ==
                 triangulation.vertices_used.size(),
               ExcInternalError());
 
        // Check whether a new level is needed. We have to check for
        // this on the highest level only
        for (const auto &cell : triangulation.active_cell_iterators_on_level(
               triangulation.levels.size() - 1))
          if (cell->refine_flag_set())
            {
              triangulation.levels.push_back(
                std::make_unique<
                  internal::TriangulationImplementation::TriaLevel>(dim));
              break;
            }
 
        // first clear user flags for quads and lines; we're going to
        // use them to flag which lines and quads need refinement
        triangulation.faces->quads.clear_user_data();
        triangulation.faces->lines.clear_user_flags();
        triangulation.faces->quads.clear_user_flags();
 
        // check how much space is needed on every level. We need not
        // check the highest level since either
        // - on the highest level no cells are flagged for refinement
        // - there are, but prepare_refinement added another empty
        //   level which then is the highest level
 
        // variables to hold the number of newly to be created
        // vertices, lines and quads. as these are stored globally,
        // declare them outside the loop over al levels. we need lines
        // and quads in pairs for refinement of old ones and lines and
        // quads, that can be stored as single ones, as they are newly
        // created in the inside of an existing cell
        unsigned int needed_vertices     = 0;
        unsigned int needed_lines_single = 0;
        unsigned int needed_quads_single = 0;
        unsigned int needed_lines_pair   = 0;
        unsigned int needed_quads_pair   = 0;
        for (int level = triangulation.levels.size() - 2; level >= 0; --level)
          {
            unsigned int new_cells = 0;
 
            for (const auto &cell :
                 triangulation.active_cell_iterators_on_level(level))
              if (cell->refine_flag_set())
                {
                  // Only support isotropic refinement
                  Assert(cell->refine_flag_set() ==
                           RefinementCase<dim>::cut_xyz,
                         ExcInternalError());
 
                  // Now count up how many new cells, faces, edges, and vertices
                  // we will need to allocate to do this refinement.
                  new_cells += cell->reference_cell().n_isotropic_children();
 
                  if (cell->reference_cell() == ReferenceCells::Hexahedron)
                    {
                      ++needed_vertices;
                      needed_lines_single += 6;
                      needed_quads_single += 12;
                    }
                  else if (cell->reference_cell() ==
                           ReferenceCells::Tetrahedron)
                    {
                      needed_lines_single += 1;
                      needed_quads_single += 8;
                    }
                  else
                    {
                      DEAL_II_ASSERT_UNREACHABLE();
                    }
 
                  // Also check whether we have to refine any of the faces and
                  // edges that bound this cell. They may of course already be
                  // refined, so we only *mark* them for refinement by setting
                  // the user flags
                  for (const auto face : cell->face_indices())
                    if (cell->face(face)->n_children() == 0)
                      cell->face(face)->set_user_flag();
                    else
                      Assert(cell->face(face)->n_children() ==
                               cell->reference_cell()
                                 .face_reference_cell(face)
                                 .n_isotropic_children(),
                             ExcInternalError());
 
                  for (const auto line : cell->line_indices())
                    if (cell->line(line)->has_children() == false)
                      cell->line(line)->set_user_flag();
                    else
                      Assert(cell->line(line)->n_children() == 2,
                             ExcInternalError());
                }
 
            const unsigned int used_cells =
              std::count(triangulation.levels[level + 1]->cells.used.begin(),
                         triangulation.levels[level + 1]->cells.used.end(),
                         true);
 
            if (triangulation.all_reference_cells_are_hyper_cube())
              reserve_space(*triangulation.levels[level + 1],
                            used_cells + new_cells,
                            3,
                            spacedim,
                            false);
            else
              reserve_space(*triangulation.levels[level + 1],
                            used_cells + new_cells,
                            3,
                            spacedim,
                            true);
 
            reserve_space(triangulation.levels[level + 1]->cells, new_cells);
          }
 
        // now count the quads and lines which were flagged for
        // refinement
        for (typename Triangulation<dim, spacedim>::quad_iterator quad =
               triangulation.begin_quad();
             quad != triangulation.end_quad();
             ++quad)
          {
            if (quad->user_flag_set() == false)
              continue;
 
            if (quad->reference_cell() == ReferenceCells::Quadrilateral)
              {
                needed_quads_pair += 4;
                needed_lines_pair += 4;
                needed_vertices += 1;
              }
            else if (quad->reference_cell() == ReferenceCells::Triangle)
              {
                needed_quads_pair += 4;
                needed_lines_single += 3;
              }
            else
              {
                DEAL_II_ASSERT_UNREACHABLE();
              }
          }
 
        for (typename Triangulation<dim, spacedim>::line_iterator line =
               triangulation.begin_line();
             line != triangulation.end_line();
             ++line)
          {
            if (line->user_flag_set() == false)
              continue;
 
            needed_lines_pair += 2;
            needed_vertices += 1;
          }
 
        reserve_space(triangulation.faces->lines,
                      needed_lines_pair,
                      needed_lines_single);
        reserve_space(*triangulation.faces,
                      needed_quads_pair,
                      needed_quads_single);
        reserve_space(triangulation.faces->quads,
                      needed_quads_pair,
                      needed_quads_single);
 
 
        // add to needed vertices how many vertices are already in use
        needed_vertices += std::count(triangulation.vertices_used.begin(),
                                      triangulation.vertices_used.end(),
                                      true);
 
        if (needed_vertices > triangulation.vertices.size())
          {
            triangulation.vertices.resize(needed_vertices, Point<spacedim>());
            triangulation.vertices_used.resize(needed_vertices, false);
          }
 
          //-----------------------------------------
          // Before we start with the actual refinement, we do some
          // sanity checks if in debug mode. especially, we try to catch
          // the notorious problem with lines being twice refined,
          // i.e. there are cells adjacent at one line ("around the
          // edge", but not at a face), with two cells differing by more
          // than one refinement level
          //
          // this check is very simple to implement here, since we have
          // all lines flagged if they shall be refined
#ifdef DEBUG
        for (const auto &cell : triangulation.active_cell_iterators())
          if (!cell->refine_flag_set())
            <