tria_accessor.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 1998 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
#ifndef dealii_tria_accessor_h
#define dealii_tria_accessor_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/bounding_box.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/point.h>
#include <deal.II/base/std_cxx26/inplace_vector.h>
#include <deal.II/base/template_constraints.h>
 
#include <deal.II/grid/cell_id.h>
#include <deal.II/grid/reference_cell.h>
#include <deal.II/grid/tria_faces.h>
#include <deal.II/grid/tria_iterator.h>
#include <deal.II/grid/tria_iterator_base.h>
#include <deal.II/grid/tria_iterator_selector.h>
#include <deal.II/grid/tria_levels.h>
 
#include <boost/container/small_vector.hpp>
 
#include <cmath>
#include <limits>
#include <set>
#include <utility>
 
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declarations
#ifndef DOXYGEN
template <int dim, int spacedim>
DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
class Triangulation;
template <typename Accessor>
class TriaRawIterator;
template <typename Accessor>
class TriaIterator;
template <typename Accessor>
class TriaActiveIterator;
 
namespace parallel
{
  template <int dim, int spacedim>
  DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
  class TriangulationBase;
}
 
template <int dim, int spacedim>
DEAL_II_CXX20_REQUIRES((concepts::is_valid_dim_spacedim<dim, spacedim>))
class DoFHandler;
template <int dim, int spacedim, bool lda>
class DoFCellAccessor;
 
 
template <int dim, int spacedim>
class Manifold;
 
template <int dim, int spacedim>
class Mapping;
#endif
 
namespace internal
{
  namespace TriangulationImplementation
  {
    class TriaObjects;
    struct Implementation;
    struct ImplementationMixedMesh;
  } // namespace TriangulationImplementation
 
  namespace TriaAccessorImplementation
  {
    struct Implementation;
 
    template <bool level_access>
    class LevelAccessor
    {};
 
    template <>
    class LevelAccessor<true>
    {
    public:
      LevelAccessor(const int level)
        : present_level(level)
      {}
 
      void
      set_level(const int level)
      {
        present_level = level;
      }
 
      int
      get_level() const
      {
        return present_level;
      }
 
    private:
      int present_level;
    };
 
    template <>
    class LevelAccessor<false>
    {
    public:
      LevelAccessor(const int level)
      {
        Assert(level == 0, ExcInternalError());
      }
 
      static void
      set_level(const int level)
      {
        Assert(level == 0, ExcInternalError());
      }
 
      static int
      get_level()
      {
        return 0;
      }
    };
  } // namespace TriaAccessorImplementation
} // namespace internal
template <int structdim, int dim, int spacedim>
class TriaAccessor;
template <int dim, int spacedim>
class TriaAccessor<0, dim, spacedim>;
template <int spacedim>
class TriaAccessor<0, 1, spacedim>;
 
namespace TriaAccessorExceptions
{
  DeclExceptionMsg(ExcCellNotUsed,
                   "The operation you are attempting can only be performed for "
                   "(cell, face, or edge) iterators that point to valid "
                   "objects. These objects need not necessarily be active, "
                   "i.e., have no children, but they need to be part of a "
                   "triangulation. (The objects pointed to by an iterator "
                   "may -- after coarsening -- also be objects that used "
                   "to be part of a triangulation, but are now no longer "
                   "used. Their memory location may have been retained "
                   "for re-use upon the next mesh refinement, but is "
                   "currently unused.)");
  DeclExceptionMsg(ExcCellNotActive,
                   "The operation you are attempting can only be performed for "
                   "(cell, face, or edge) iterators that point to 'active' "
                   "objects. 'Active' objects are those that do not have "
                   "children (in the case of cells), or that are part of "
                   "an active cell (in the case of faces or edges). However, "
                   "the object on which you are trying the current "
                   "operation is not 'active' in this sense.");
  DeclExceptionMsg(ExcCellHasNoChildren,
                   "The operation you are attempting can only be performed for "
                   "(cell, face, or edge) iterators that have children, "
                   "but the object on which you are trying the current "
                   "operation does not have any.");
  DeclExceptionMsg(ExcCellHasNoParent,
                   "The operation you are attempting can only be performed for "
                   "(cell, face, or edge) iterators that have a parent object, "
                   "but the object on which you are trying the current "
                   "operation does not have one -- i.e., it is on the "
                   "coarsest level of the triangulation.");
  DeclException1(ExcCantSetChildren,
                 int,
                 << "You can only set the child index if the cell does not "
                 << "currently have children registered; or you can clear it. "
                 << "The given index was " << arg1
                 << " (-1 means: clear children).");
  template <typename AccessorType>
  DeclException1(ExcDereferenceInvalidObject,
                 AccessorType,
                 << "You tried to dereference an iterator for which this "
                 << "is not possible. More information on this iterator: "
                 << "index=" << arg1.index() << ", state="
                 << (arg1.state() == IteratorState::valid ?
                       "valid" :
                       (arg1.state() == IteratorState::past_the_end ?
                          "past_the_end" :
                          "invalid")));
  DeclExceptionMsg(ExcCantCompareIterators,
                   "Iterators can only be compared if they point to the same "
                   "triangulation, or if neither of them are associated "
                   "with a triangulation.");
  // TODO: Write documentation!
  DeclException0(ExcNeighborIsCoarser);
  // TODO: Write documentation!
  DeclException0(ExcNeighborIsNotCoarser);
  DeclException0(ExcFacesHaveNoLevel);
  DeclException0(ExcNoPeriodicNeighbor);
  DeclException1(
    ExcSetOnlyEvenChildren,
    int,
    << "You can only set the child index of an even numbered child."
    << "The number of the child given was " << arg1 << '.');
} // namespace TriaAccessorExceptions
 
 
template <int structdim, int dim, int spacedim = dim>
class TriaAccessorBase
  : private internal::TriaAccessorImplementation::LevelAccessor<structdim ==
                                                                dim>
{
public:
  static constexpr unsigned int space_dimension = spacedim;
 
  static constexpr unsigned int dimension = dim;
 
  static const unsigned int structure_dimension = structdim;
 
  void
  operator=(const TriaAccessorBase *) = delete;
 
protected:
  using AccessorData = void;
 
  TriaAccessorBase(const Triangulation<dim, spacedim> *parent = nullptr,
                   const int                           level  = -1,
                   const int                           index  = -1,
                   const AccessorData                       * = nullptr);
 
  TriaAccessorBase(const TriaAccessorBase &);
 
  void
  copy_from(const TriaAccessorBase &);
 
  TriaAccessorBase &
  operator=(const TriaAccessorBase &);
 
  bool
  operator<(const TriaAccessorBase &other) const;
 
protected:
  bool
  operator==(const TriaAccessorBase &) const;
 
  bool
  operator!=(const TriaAccessorBase &) const;
 
  void
  operator++();
 
  void
  operator--();
  ::internal::TriangulationImplementation::TriaObjects &
  objects() const;
 
public:
  using LocalData = void *;
 
  int
  level() const;
 
  int
  index() const;
 
  IteratorState::IteratorStates
  state() const;
 
  const Triangulation<dim, spacedim> &
  get_triangulation() const;
 
protected:
  int present_index;
 
  const Triangulation<dim, spacedim> *tria;
 
private:
  template <typename Accessor>
  friend class TriaRawIterator;
  template <typename Accessor>
  friend class TriaIterator;
  template <typename Accessor>
  friend class TriaActiveIterator;
};
 
 
 
template <int structdim, int dim, int spacedim = dim>
class InvalidAccessor
{
public:
  static constexpr unsigned int space_dimension = spacedim;
 
  static constexpr unsigned int dimension = dim;
 
  static const unsigned int structure_dimension = structdim;
 
  using AccessorData = void;
 
  InvalidAccessor(const void         *parent     = nullptr,
                  const int           level      = -1,
                  const int           index      = -1,
                  const AccessorData *local_data = nullptr);
 
  InvalidAccessor(const InvalidAccessor &);
 
  template <typename OtherAccessor>
  InvalidAccessor(const OtherAccessor &);
 
  void
  copy_from(const InvalidAccessor &);
 
  bool
  operator==(const InvalidAccessor &) const;
  bool
  operator!=(const InvalidAccessor &) const;
 
  void
  operator++() const;
  void
  operator--() const;
 
  static IteratorState::IteratorStates
  state();
 
 
  static int
  level();
 
  static int
  index();
 
  bool
  used() const;
 
  bool
  has_children() const;
 
  types::manifold_id
  manifold_id() const;
 
  unsigned int
  user_index() const;
 
  void
  set_user_index(const unsigned int p) const;
 
  void
  set_manifold_id(const types::manifold_id) const;
 
  Point<spacedim> &
  vertex(const unsigned int i) const;
 
  void *
  line(const unsigned int i) const;
 
  void *
  quad(const unsigned int i) const;
};
 
 
 
template <int structdim, int dim, int spacedim>
class TriaAccessor : public TriaAccessorBase<structdim, dim, spacedim>
{
public:
  using AccessorData =
    typename TriaAccessorBase<structdim, dim, spacedim>::AccessorData;
 
  TriaAccessor(const Triangulation<dim, spacedim> *parent     = nullptr,
               const int                           level      = -1,
               const int                           index      = -1,
               const AccessorData                 *local_data = nullptr);
 
  TriaAccessor(const TriaAccessor &) = default;
 
  TriaAccessor(TriaAccessor &&) = default; // NOLINT
 
  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor(const InvalidAccessor<structdim2, dim2, spacedim2> &);
 
  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor(const TriaAccessor<structdim2, dim2, spacedim2> &);
 
  TriaAccessor &
  operator=(const TriaAccessor &) = delete;
 
  TriaAccessor &
  operator=(TriaAccessor &&) = default; // NOLINT
 
  ~TriaAccessor() = default;
 
  bool
  used() const;
 
  TriaIterator<TriaAccessor<0, dim, spacedim>>
  vertex_iterator(const unsigned int i) const;
 
  unsigned int
  vertex_index(const unsigned int i) const;
 
  Point<spacedim> &
  vertex(const unsigned int i) const;
 
  typename ::internal::TriangulationImplementation::
    Iterators<dim, spacedim>::line_iterator
    line(const unsigned int i) const;
 
  unsigned int
  line_index(const unsigned int i) const;
 
  typename ::internal::TriangulationImplementation::
    Iterators<dim, spacedim>::quad_iterator
    quad(const unsigned int i) const;
 
  unsigned int
  quad_index(const unsigned int i) const;
  types::geometric_orientation
  combined_face_orientation(const unsigned int face) const;
 
  bool
  face_orientation(const unsigned int face) const;
 
  bool
  face_flip(const unsigned int face) const;
 
  bool
  face_rotation(const unsigned int face) const;
 
  types::geometric_orientation
  line_orientation(const unsigned int line) const;
  bool
  has_children() const;
 
  unsigned int
  n_children() const;
 
  unsigned int
  n_active_descendants() const;
 
  unsigned int
  max_refinement_depth() const;
 
  TriaIterator<TriaAccessor<structdim, dim, spacedim>>
  child(const unsigned int i) const;
 
  unsigned int
  child_iterator_to_index(
    const TriaIterator<TriaAccessor<structdim, dim, spacedim>> &child) const;
 
  TriaIterator<TriaAccessor<structdim, dim, spacedim>>
  isotropic_child(const unsigned int i) const;
 
  RefinementCase<structdim>
  refinement_case() const;
 
  int
  child_index(const unsigned int i) const;
 
  int
  isotropic_child_index(const unsigned int i) const;
  types::boundary_id
  boundary_id() const;
 
  void
  set_boundary_id(const types::boundary_id) const;
 
  void
  set_all_boundary_ids(const types::boundary_id) const;
 
  bool
  at_boundary() const;
 
  const Manifold<dim, spacedim> &
  get_manifold() const;
 
  types::manifold_id
  manifold_id() const;
 
  void
  set_manifold_id(const types::manifold_id) const;
 
  void
  set_all_manifold_ids(const types::manifold_id) const;
 
  bool
  user_flag_set() const;
 
  void
  set_user_flag() const;
 
  void
  clear_user_flag() const;
 
  void
  recursively_set_user_flag() const;
 
  void
  recursively_clear_user_flag() const;
 
  void
  clear_user_data() const;
 
  void
  set_user_pointer(void *p) const;
 
  void
  clear_user_pointer() const;
 
  void *
  user_pointer() const;
 
  void
  recursively_set_user_pointer(void *p) const;
 
  void
  recursively_clear_user_pointer() const;
 
  void
  set_user_index(const unsigned int p) const;
 
  void
  clear_user_index() const;
 
  unsigned int
  user_index() const;
 
  void
  recursively_set_user_index(const unsigned int p) const;
 
  void
  recursively_clear_user_index() const;
  double
  diameter() const;
 
  std::pair<Point<spacedim>, double>
  enclosing_ball() const;
 
  BoundingBox<spacedim>
  bounding_box() const;
 
  double
  extent_in_direction(const unsigned int axis) const;
 
  double
  minimum_vertex_distance() const;
 
  Point<spacedim>
  intermediate_point(const Point<structdim> &coordinates) const;
 
  Point<structdim>
  real_to_unit_cell_affine_approximation(const Point<spacedim> &point) const;
 
  Point<spacedim>
  center(const bool respect_manifold             = false,
         const bool interpolate_from_surrounding = false) const;
 
  Point<spacedim>
  barycenter() const;
 
  double
  measure() const;
 
  bool
  is_translation_of(
    const TriaIterator<TriaAccessor<structdim, dim, spacedim>> &o) const;
 
  ReferenceCell<structdim>
  reference_cell() const;
 
  unsigned int
  n_vertices() const;
 
  unsigned int
  n_lines() const;
 
  unsigned int
  n_faces() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  vertex_indices() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  line_indices() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  face_indices() const;
 
private:
  void
  set_boundary_id_internal(const types::boundary_id id) const;
 
  void
  set_bounding_object_indices(
    const std::initializer_list<int> &new_indices) const;
 
  void
  set_bounding_object_indices(
    const std::initializer_list<unsigned int> &new_indices) const;
 
  void
  set_line_orientation(const unsigned int                 line,
                       const types::geometric_orientation orientation) const;
 
  void
  set_combined_face_orientation(
    const unsigned int                 face_no,
    const types::geometric_orientation combined_orientation) const;
 
  void
  set_used_flag() const;
 
  void
  clear_used_flag() const;
 
  void
  set_refinement_case(const RefinementCase<structdim> &ref_case) const;
 
  void
  clear_refinement_case() const;
 
  void
  set_children(const unsigned int i, const int index) const;
 
  void
  clear_children() const;
 
private:
  friend class Triangulation<dim, spacedim>;
 
  friend struct ::internal::TriangulationImplementation::Implementation;
  friend struct ::internal::TriangulationImplementation::
    ImplementationMixedMesh;
  friend struct ::internal::TriaAccessorImplementation::Implementation;
};
 
 
 
template <int dim, int spacedim>
class TriaAccessor<0, dim, spacedim>
{
public:
  static constexpr unsigned int space_dimension = spacedim;
 
  static constexpr unsigned int dimension = dim;
 
  static const unsigned int structure_dimension = 0;
 
  using AccessorData = void;
 
  TriaAccessor(const Triangulation<dim, spacedim> *tria,
               const unsigned int                  vertex_index);
 
  TriaAccessor(const Triangulation<dim, spacedim> *tria  = nullptr,
               const int                           level = 0,
               const int                           index = 0,
               const AccessorData                      * = nullptr);
 
  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor(const TriaAccessor<structdim2, dim2, spacedim2> &);
 
  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor(const InvalidAccessor<structdim2, dim2, spacedim2> &);
 
  IteratorState::IteratorStates
  state() const;
 
  static int
  level();
 
  int
  index() const;
 
  const Triangulation<dim, spacedim> &
  get_triangulation() const;
 
  void
  operator++();
 
  void
  operator--();
  bool
  operator==(const TriaAccessor &) const;
 
  bool
  operator!=(const TriaAccessor &) const;
 
  unsigned int
  vertex_index(const unsigned int i = 0) const;
 
  Point<spacedim> &
  vertex(const unsigned int i = 0) const;
 
  typename ::internal::TriangulationImplementation::
    Iterators<dim, spacedim>::line_iterator static line(const unsigned int);
 
  static unsigned int
  line_index(const unsigned int i);
 
  static typename ::internal::TriangulationImplementation::
    Iterators<dim, spacedim>::quad_iterator
    quad(const unsigned int i);
 
  static unsigned int
  quad_index(const unsigned int i);
 
  double
  diameter() const;
 
  double
  extent_in_direction(const unsigned int axis) const;
 
  Point<spacedim>
  center(const bool respect_manifold             = false,
         const bool interpolate_from_surrounding = false) const;
 
  double
  measure() const;
  static types::geometric_orientation
  combined_face_orientation(const unsigned int face);
 
  static bool
  face_orientation(const unsigned int face);
 
  static bool
  face_flip(const unsigned int face);
 
  static bool
  face_rotation(const unsigned int face);
 
  static types::geometric_orientation
  line_orientation(const unsigned int line);
 
  static bool
  has_children();
 
  static unsigned int
  n_children();
 
  static unsigned int
  n_active_descendants();
 
  static unsigned int
  max_refinement_depth();
 
  static unsigned int
  child_iterator_to_index(const TriaIterator<TriaAccessor<0, dim, spacedim>> &);
 
  static TriaIterator<TriaAccessor<0, dim, spacedim>>
  child(const unsigned int);
 
  static TriaIterator<TriaAccessor<0, dim, spacedim>>
  isotropic_child(const unsigned int);
 
  static RefinementCase<0>
  refinement_case();
 
  static int
  child_index(const unsigned int i);
 
  static int
  isotropic_child_index(const unsigned int i);
  bool
  used() const;
 
protected:
  void
  copy_from(const TriaAccessor &);
 
  bool
  operator<(const TriaAccessor &other) const;
 
  const Triangulation<dim, spacedim> *tria;
 
  unsigned int global_vertex_index;
 
private:
  template <typename Accessor>
  friend class TriaRawIterator;
  template <typename Accessor>
  friend class TriaIterator;
  template <typename Accessor>
  friend class TriaActiveIterator;
};
 
 
 
template <int spacedim>
class TriaAccessor<0, 1, spacedim>
{
public:
  static constexpr unsigned int space_dimension = spacedim;
 
  static constexpr unsigned int dimension = 1;
 
  static const unsigned int structure_dimension = 0;
 
  using AccessorData = void;
 
  enum VertexKind
  {
    left_vertex,
    interior_vertex,
    right_vertex
  };
 
  TriaAccessor(const Triangulation<1, spacedim> *tria,
               const VertexKind                  vertex_kind,
               const unsigned int                vertex_index);
 
  TriaAccessor(const Triangulation<1, spacedim> *tria = nullptr,
               const int                              = 0,
               const int                              = 0,
               const AccessorData                   * = nullptr);
 
  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor(const TriaAccessor<structdim2, dim2, spacedim2> &);
 
  template <int structdim2, int dim2, int spacedim2>
  TriaAccessor(const InvalidAccessor<structdim2, dim2, spacedim2> &);
 
  void
  copy_from(const TriaAccessor &);
 
  void
  copy_from(const TriaAccessorBase<0, 1, spacedim> &);
 
  static IteratorState::IteratorStates
  state();
 
  static int
  level();
 
  int
  index() const;
 
  const Triangulation<1, spacedim> &
  get_triangulation() const;
 
  void
  operator++() const;
 
  void
  operator--() const;
  bool
  operator==(const TriaAccessor &) const;
 
  bool
  operator!=(const TriaAccessor &) const;
 
  bool
  operator<(const TriaAccessor &other) const;
 
  unsigned int
  vertex_index(const unsigned int i = 0) const;
 
  Point<spacedim> &
  vertex(const unsigned int i = 0) const;
 
  Point<spacedim>
  center() const;
 
  typename ::internal::TriangulationImplementation::
    Iterators<1, spacedim>::line_iterator static line(const unsigned int);
 
  static unsigned int
  line_index(const unsigned int i);
 
  static typename ::internal::TriangulationImplementation::
    Iterators<1, spacedim>::quad_iterator
    quad(const unsigned int i);
 
  static unsigned int
  quad_index(const unsigned int i);
 
  static double
  measure();
 
  bool
  at_boundary() const;
 
  types::boundary_id
  boundary_id() const;
 
  const Manifold<1, spacedim> &
  get_manifold() const;
 
  types::manifold_id
  manifold_id() const;
 
 
  bool
  user_flag_set() const;
 
  void
  set_user_flag() const;
 
  void
  clear_user_flag() const;
 
  void
  recursively_set_user_flag() const;
 
  void
  recursively_clear_user_flag() const;
 
  void
  clear_user_data() const;
 
  void
  set_user_pointer(void *p) const;
 
  void
  clear_user_pointer() const;
 
  void *
  user_pointer() const;
 
  void
  recursively_set_user_pointer(void *p) const;
 
  void
  recursively_clear_user_pointer() const;
 
  void
  set_user_index(const unsigned int p) const;
 
  void
  clear_user_index() const;
 
  unsigned int
  user_index() const;
 
  void
  recursively_set_user_index(const unsigned int p) const;
 
  void
  recursively_clear_user_index() const;
  static types::geometric_orientation
  combined_face_orientation(const unsigned int face);
 
  static bool
  face_orientation(const unsigned int face);
 
  static bool
  face_flip(const unsigned int face);
 
  static bool
  face_rotation(const unsigned int face);
 
  static types::geometric_orientation
  line_orientation(const unsigned int line);
 
  static bool
  has_children();
 
  static unsigned int
  n_children();
 
  static unsigned int
  n_active_descendants();
 
  static unsigned int
  max_refinement_depth();
 
  static unsigned int
  child_iterator_to_index(const TriaIterator<TriaAccessor<0, 1, spacedim>> &);
 
  static TriaIterator<TriaAccessor<0, 1, spacedim>>
  child(const unsigned int);
 
  static TriaIterator<TriaAccessor<0, 1, spacedim>>
  isotropic_child(const unsigned int);
 
  static RefinementCase<0>
  refinement_case();
 
  static int
  child_index(const unsigned int i);
 
  static int
  isotropic_child_index(const unsigned int i);
  void
  set_boundary_id(const types::boundary_id) const;
 
  void
  set_manifold_id(const types::manifold_id);
 
  void
  set_all_boundary_ids(const types::boundary_id) const;
 
  void
  set_all_manifold_ids(const types::manifold_id);
  bool
  used() const;
 
  ReferenceCell<0>
  reference_cell() const;
 
  unsigned int
  n_vertices() const;
 
  unsigned int
  n_lines() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  vertex_indices() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  line_indices() const;
 
protected:
  const Triangulation<1, spacedim> *tria;
 
  VertexKind vertex_kind;
 
  unsigned int global_vertex_index;
};
 
 
 
template <int dim, int spacedim = dim>
class CellAccessor : public TriaAccessor<dim, dim, spacedim>
{
public:
  using AccessorData = typename TriaAccessor<dim, dim, spacedim>::AccessorData;
 
  using Container = Triangulation<dim, spacedim>;
 
  CellAccessor(const Triangulation<dim, spacedim> *parent     = nullptr,
               const int                           level      = -1,
               const int                           index      = -1,
               const AccessorData                 *local_data = nullptr);
 
  CellAccessor(const TriaAccessor<dim, dim, spacedim> &cell_accessor);
 
  template <int structdim2, int dim2, int spacedim2>
  CellAccessor(const InvalidAccessor<structdim2, dim2, spacedim2> &);
 
  template <int structdim2, int dim2, int spacedim2>
  CellAccessor(const TriaAccessor<structdim2, dim2, spacedim2> &);
 
  CellAccessor(const CellAccessor<dim, spacedim> &) = default;
 
  // NOLINTNEXTLINE OSX does not compile with noexcept
  CellAccessor(CellAccessor<dim, spacedim> &&) = default;
 
  ~CellAccessor() = default;
 
  CellAccessor<dim, spacedim> &
  operator=(const CellAccessor<dim, spacedim> &) = delete;
 
  // NOLINTNEXTLINE OSX does not compile with noexcept
  CellAccessor<dim, spacedim> &
  operator=(CellAccessor<dim, spacedim> &&) = default; // NOLINT
 
  TriaActiveIterator<DoFCellAccessor<dim, spacedim, false>>
  as_dof_handler_iterator(const DoFHandler<dim, spacedim> &dof_handler) const;
 
  TriaIterator<DoFCellAccessor<dim, spacedim, true>>
  as_dof_handler_level_iterator(
    const DoFHandler<dim, spacedim> &dof_handler) const;
 
 
  TriaIterator<CellAccessor<dim, spacedim>>
  child(const unsigned int i) const;
 
  boost::container::small_vector<TriaIterator<CellAccessor<dim, spacedim>>,
                                 ReferenceCells::max_n_children<dim>()>
  child_iterators() const;
 
  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
  face(const unsigned int i) const;
 
  unsigned int
  face_iterator_to_index(
    const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> &face) const;
 
  boost::container::small_vector<
    TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
#ifndef _MSC_VER // MSVC prior to 2022 cannot use a constexpr function this way
    ReferenceCells::max_n_faces<dim>()
#else
    GeometryInfo<dim>::faces_per_cell
#endif
    >
  face_iterators() const;
 
  unsigned int
  face_index(const unsigned int i) const;
 
  TriaIterator<CellAccessor<dim, spacedim>>
  neighbor_child_on_subface(const unsigned int face_no,
                            const unsigned int subface_no) const;
 
  TriaIterator<CellAccessor<dim, spacedim>>
  neighbor(const unsigned int face_no) const;
 
  int
  neighbor_index(const unsigned int face_no) const;
 
  int
  neighbor_level(const unsigned int face_no) const;
 
  unsigned int
  neighbor_of_neighbor(const unsigned int face_no) const;
 
  bool
  neighbor_is_coarser(const unsigned int face_no) const;
 
  std::pair<unsigned int, unsigned int>
  neighbor_of_coarser_neighbor(const unsigned int neighbor) const;
 
  unsigned int
  neighbor_face_no(const unsigned int neighbor) const;
 
  static bool
  is_level_cell();
 
  bool
  has_periodic_neighbor(const unsigned int i) const;
 
  TriaIterator<CellAccessor<dim, spacedim>>
  periodic_neighbor(const unsigned int i) const;
 
  TriaIterator<CellAccessor<dim, spacedim>>
  neighbor_or_periodic_neighbor(const unsigned int i) const;
 
  TriaIterator<CellAccessor<dim, spacedim>>
  periodic_neighbor_child_on_subface(const unsigned int face_no,
                                     const unsigned int subface_no) const;
 
  std::pair<unsigned int, unsigned int>
  periodic_neighbor_of_coarser_periodic_neighbor(const unsigned face_no) const;
 
  int
  periodic_neighbor_index(const unsigned int i) const;
 
  int
  periodic_neighbor_level(const unsigned int i) const;
 
  unsigned int
  periodic_neighbor_of_periodic_neighbor(const unsigned int i) const;
 
  unsigned int
  periodic_neighbor_face_no(const unsigned int i) const;
 
  bool
  periodic_neighbor_is_coarser(const unsigned int i) const;
 
  bool
  at_boundary(const unsigned int i) const;
 
  bool
  at_boundary() const;
 
  bool
  has_boundary_lines() const;
  RefinementCase<dim>
  refine_flag_set() const;
 
  void
  set_refine_flag(const RefinementCase<dim> ref_case =
                    RefinementCase<dim>::isotropic_refinement) const;
 
  void
  clear_refine_flag() const;
 
  std::uint8_t
  refine_choice() const;
 
  void
  set_refine_choice(const std::uint8_t refinement_choice = static_cast<char>(
                      IsotropicRefinementChoice::isotropic_refinement)) const;
 
  void
  clear_refine_choice() const;
 
  bool
  flag_for_face_refinement(
    const unsigned int             face_no,
    const RefinementCase<dim - 1> &face_refinement_case =
      RefinementCase<dim - 1>::isotropic_refinement) const;
 
  bool
  flag_for_line_refinement(const unsigned int line_no) const;
 
  internal::SubfaceCase<dim>
  subface_case(const unsigned int face_no) const;
 
  bool
  coarsen_flag_set() const;
 
  void
  set_coarsen_flag() const;
 
  void
  clear_coarsen_flag() const;
  types::material_id
  material_id() const;
 
  void
  set_material_id(const types::material_id new_material_id) const;
 
  void
  recursively_set_material_id(const types::material_id new_material_id) const;
  types::subdomain_id
  subdomain_id() const;
 
  void
  set_subdomain_id(const types::subdomain_id new_subdomain_id) const;
 
  types::subdomain_id
  level_subdomain_id() const;
 
  void
  set_level_subdomain_id(
    const types::subdomain_id new_level_subdomain_id) const;
 
 
  void
  recursively_set_subdomain_id(
    const types::subdomain_id new_subdomain_id) const;
  types::global_cell_index
  global_active_cell_index() const;
 
  types::global_cell_index
  global_level_cell_index() const;
 
  unsigned int
  active_cell_index() const;
 
  int
  parent_index() const;
 
  TriaIterator<CellAccessor<dim, spacedim>>
  parent() const;
 
  bool
  direction_flag() const;
 
  bool
  is_active() const;
 
  bool
  is_locally_owned() const;
 
  bool
  is_locally_owned_on_level() const;
 
  bool
  is_ghost() const;
 
  bool
  is_ghost_on_level() const;
 
  bool
  is_artificial() const;
 
  bool
  is_artificial_on_level() const;
 
  bool
  point_inside(const Point<spacedim> &p) const;
 
  DEAL_II_DEPRECATED
  std::set<TriaActiveIterator<CellAccessor<dim, spacedim>>>
  get_cells_adjacent_to_line(const unsigned int i) const;
 
  CellId
  id() const;
 
  using TriaAccessor<dim, dim, spacedim>::diameter;
 
  double
  diameter(const Mapping<dim, spacedim> &mapping) const;
 
  DeclException0(ExcRefineCellNotActive);
  DeclException0(ExcCellFlaggedForRefinement);
  DeclException0(ExcCellFlaggedForCoarsening);
 
protected:
  void
  set_neighbor(const unsigned int                               face_no,
               const TriaIterator<CellAccessor<dim, spacedim>> &neighbor) const;
 
  unsigned int
  neighbor_of_neighbor_internal(const unsigned int neighbor) const;
 
  template <int dim_, int spacedim_>
  bool
  point_inside_codim(const Point<spacedim_> &p) const;
 
 
 
private:
  void
  set_active_cell_index(const unsigned int active_cell_index) const;
 
  void
  set_global_active_cell_index(const types::global_cell_index index) const;
 
  void
  set_global_level_cell_index(const types::global_cell_index index) const;
 
  void
  set_parent(const unsigned int parent_index);
 
  void
  set_direction_flag(const bool new_direction_flag) const;
 
  friend class Triangulation<dim, spacedim>;
 
  friend class parallel::TriangulationBase<dim, spacedim>;
 
  friend struct ::internal::TriangulationImplementation::Implementation;
  friend struct ::internal::TriangulationImplementation::
    ImplementationMixedMesh;
};
 
 
 
/* ----- declaration of explicit specializations and general templates ----- */
 
 
template <int structdim, int dim, int spacedim>
template <typename OtherAccessor>
InvalidAccessor<structdim, dim, spacedim>::InvalidAccessor(
  const OtherAccessor &)
{
  Assert(false,
         ExcMessage("You are attempting an illegal conversion between "
                    "iterator/accessor types. The constructor you call "
                    "only exists to make certain template constructs "
                    "easier to write as dimension independent code but "
                    "the conversion is not valid in the current context."));
}
 
 
 
template <int structdim, int dim, int spacedim>
template <int structdim2, int dim2, int spacedim2>
TriaAccessor<structdim, dim, spacedim>::TriaAccessor(
  const InvalidAccessor<structdim2, dim2, spacedim2> &)
{
  Assert(false,
         ExcMessage("You are attempting an illegal conversion between "
                    "iterator/accessor types. The constructor you call "
                    "only exists to make certain template constructs "
                    "easier to write as dimension independent code but "
                    "the conversion is not valid in the current context."));
}
 
 
 
template <int dim, int spacedim>
template <int structdim2, int dim2, int spacedim2>
CellAccessor<dim, spacedim>::CellAccessor(
  const InvalidAccessor<structdim2, dim2, spacedim2> &)
{
  Assert(false,
         ExcMessage("You are attempting an illegal conversion between "
                    "iterator/accessor types. The constructor you call "
                    "only exists to make certain template constructs "
                    "easier to write as dimension independent code but "
                    "the conversion is not valid in the current context."));
}
 
 
 
template <int structdim, int dim, int spacedim>
template <int structdim2, int dim2, int spacedim2>
TriaAccessor<structdim, dim, spacedim>::TriaAccessor(
  const TriaAccessor<structdim2, dim2, spacedim2> &)
{
  Assert(false,
         ExcMessage("You are attempting an illegal conversion between "
                    "iterator/accessor types. The constructor you call "
                    "only exists to make certain template constructs "
                    "easier to write as dimension independent code but "
                    "the conversion is not valid in the current context."));
}
 
 
 
template <int dim, int spacedim>
template <int structdim2, int dim2, int spacedim2>
CellAccessor<dim, spacedim>::CellAccessor(
  const TriaAccessor<structdim2, dim2, spacedim2> &)
{
  Assert(false,
         ExcMessage("You are attempting an illegal conversion between "
                    "iterator/accessor types. The constructor you call "
                    "only exists to make certain template constructs "
                    "easier to write as dimension independent code but "
                    "the conversion is not valid in the current context."));
}
 
 
#ifndef DOXYGEN
 
template <>
bool
CellAccessor<1, 1>::point_inside(const Point<1> &) const;
template <>
bool
CellAccessor<2, 2>::point_inside(const Point<2> &) const;
template <>
bool
CellAccessor<3, 3>::point_inside(const Point<3> &) const;
template <>
bool
CellAccessor<1, 2>::point_inside(const Point<2> &) const;
template <>
bool
CellAccessor<1, 3>::point_inside(const Point<3> &) const;
template <>
bool
CellAccessor<2, 3>::point_inside(const Point<3> &) const;
// -------------------------------------------------------------------
 
template <>
void
TriaAccessor<3, 3, 3>::set_all_manifold_ids(const types::manifold_id) const;
 
 
 
namespace internal
{
  namespace TriaAccessorImplementation
  {
    template <int dim, int spacedim>
    inline double
    diameter(const ArrayView<const Point<spacedim>> &vertices)
    {
      const ReferenceCell<dim> reference_cell =
        ReferenceCells::n_vertices_to_reference_cell<dim>(vertices.size());
 
      if constexpr (dim == 1)
        {
          Assert(reference_cell == ReferenceCells::Line, ExcInternalError());
          // Return the distance between the two vertices
          return (vertices[1] - vertices[0]).norm();
        }
      else if constexpr (dim == 2)
        {
          if (reference_cell == ReferenceCells::Triangle)
            // Return the longest of the three edges
            return std::max({(vertices[1] - vertices[0]).norm(),
                             (vertices[2] - vertices[1]).norm(),
                             (vertices[2] - vertices[0]).norm()});
          else if (reference_cell == ReferenceCells::Quadrilateral)
            // Return the longer one of the two diagonals of the quadrilateral
            return std::max({(vertices[3] - vertices[0]).norm(),
                             (vertices[2] - vertices[1]).norm()});
        }
      else if constexpr (dim == 3)
        {
          if (reference_cell == ReferenceCells::Tetrahedron)
            // Return the longest of the six edges of the tetrahedron
            return std::max({(vertices[1] - vertices[0]).norm(),
                             (vertices[2] - vertices[0]).norm(),
                             (vertices[2] - vertices[1]).norm(),
                             (vertices[3] - vertices[0]).norm(),
                             (vertices[3] - vertices[1]).norm(),
                             (vertices[3] - vertices[2]).norm()});
          else if (reference_cell == ReferenceCells::Pyramid)
            // Return ...
            return std::max({// the longest diagonal of the quadrilateral base
                             // of the pyramid or ...
                             (vertices[3] - vertices[0]).norm(),
                             (vertices[2] - vertices[1]).norm(),
                             // the longest edge connected with the apex of the
                             // pyramid
                             (vertices[4] - vertices[0]).norm(),
                             (vertices[4] - vertices[1]).norm(),
                             (vertices[4] - vertices[2]).norm(),
                             (vertices[4] - vertices[3]).norm()});
          else if (reference_cell == ReferenceCells::Wedge)
            // Return ...
            return std::max({// the longest of the 2*3=6 diagonals of the three
                             // quadrilateral sides of the wedge or ...
                             (vertices[4] - vertices[0]).norm(),
                             (vertices[3] - vertices[1]).norm(),
                             (vertices[5] - vertices[1]).norm(),
                             (vertices[4] - vertices[2]).norm(),
                             (vertices[5] - vertices[0]).norm(),
                             (vertices[3] - vertices[2]).norm(),
                             // the longest of the 3*2=6 edges of the two
                             // triangular faces of the wedge
                             (vertices[1] - vertices[0]).norm(),
                             (vertices[2] - vertices[1]).norm(),
                             (vertices[2] - vertices[0]).norm(),
                             (vertices[4] - vertices[3]).norm(),
                             (vertices[5] - vertices[4]).norm(),
                             (vertices[5] - vertices[3]).norm()});
          else if (reference_cell == ReferenceCells::Hexahedron)
            // Return the longest of the four diagonals of the hexahedron
            return std::max({(vertices[7] - vertices[0]).norm(),
                             (vertices[6] - vertices[1]).norm(),
                             (vertices[2] - vertices[5]).norm(),
                             (vertices[3] - vertices[4]).norm()});
        }
 
      DEAL_II_NOT_IMPLEMENTED();
      return -1e10;
    }
  } // namespace TriaAccessorImplementation
} // namespace internal
 
 
/*--------------------- Functions: TriaAccessorBase -------------------------*/
 
template <int structdim, int dim, int spacedim>
inline TriaAccessorBase<structdim, dim, spacedim>::TriaAccessorBase(
  const Triangulation<dim, spacedim> *tria,
  const int                           level,
  const int                           index,
  const AccessorData *)
  : internal::TriaAccessorImplementation::LevelAccessor<dim == structdim>(
      (structdim == dim) ? level : 0)
  , present_index(index)
  , tria(tria)
{
  // non-cells have no level, so a 0 should have been passed, or a -1 for an
  // end-iterator, or -2 for an invalid (default constructed) iterator
  if (structdim != dim)
    Assert((level == 0) || (level == -1) || (level == -2), ExcInternalError());
}
 
 
template <int structdim, int dim, int spacedim>
inline TriaAccessorBase<structdim, dim, spacedim>::TriaAccessorBase(
  const TriaAccessorBase<structdim, dim, spacedim> &a)
  : internal::TriaAccessorImplementation::LevelAccessor<dim == structdim>(
      a.level())
  , present_index(a.present_index)
  , tria(a.tria)
{}
 
 
template <int structdim, int dim, int spacedim>
inline void
TriaAccessorBase<structdim, dim, spacedim>::copy_from(
  const TriaAccessorBase<structdim, dim, spacedim> &a)
{
  this->set_level(a.level());
  present_index = a.present_index;
  tria          = a.tria;
 
  if (structdim != dim)
    Assert(level() == 0, ExcInternalError());
}
 
 
 
template <int structdim, int dim, int spacedim>
inline TriaAccessorBase<structdim, dim, spacedim> &
TriaAccessorBase<structdim, dim, spacedim>::operator=(
  const TriaAccessorBase<structdim, dim, spacedim> &a)
{
  this->set_level(a.level());
  present_index = a.present_index;
  tria          = a.tria;
 
  if (structdim != dim)
    Assert(level() == 0, ExcInternalError());
  return *this;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline bool
TriaAccessorBase<structdim, dim, spacedim>::operator==(
  const TriaAccessorBase<structdim, dim, spacedim> &a) const
{
  Assert(tria == a.tria || tria == nullptr || a.tria == nullptr,
         TriaAccessorExceptions::ExcCantCompareIterators());
  return ((tria == a.tria) && (level() == a.level()) &&
          (present_index == a.present_index));
}
 
 
 
template <int structdim, int dim, int spacedim>
inline bool
TriaAccessorBase<structdim, dim, spacedim>::operator!=(
  const TriaAccessorBase<structdim, dim, spacedim> &a) const
{
  Assert(tria == a.tria || tria == nullptr || a.tria == nullptr,
         TriaAccessorExceptions::ExcCantCompareIterators());
  return ((tria != a.tria) || (level() != a.level()) ||
          (present_index != a.present_index));
}
 
 
 
template <int structdim, int dim, int spacedim>
inline bool
TriaAccessorBase<structdim, dim, spacedim>::operator<(
  const TriaAccessorBase<structdim, dim, spacedim> &other) const
{
  Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators());
 
  if (level() != other.level())
    return (level() < other.level());
 
  return (present_index < other.present_index);
}
 
 
 
template <int structdim, int dim, int spacedim>
inline int
TriaAccessorBase<structdim, dim, spacedim>::level() const
{
  // This is always zero or invalid
  // if the object is not a cell
  return this->get_level();
}
 
 
 
template <int structdim, int dim, int spacedim>
inline int
TriaAccessorBase<structdim, dim, spacedim>::index() const
{
  return present_index;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline IteratorState::IteratorStates
TriaAccessorBase<structdim, dim, spacedim>::state() const
{
  if ((level() >= 0) && (present_index >= 0))
    return IteratorState::valid;
  else if (present_index == -1)
    return IteratorState::past_the_end;
  else
    return IteratorState::invalid;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline const Triangulation<dim, spacedim> &
TriaAccessorBase<structdim, dim, spacedim>::get_triangulation() const
{
  return *tria;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline void
TriaAccessorBase<structdim, dim, spacedim>::operator++()
{
  // this iterator is used for
  // objects without level
  ++this->present_index;
 
  if (structdim != dim)
    {
      // is index still in the range of
      // the vector? (note that we don't
      // have to set the level, since
      // dim!=1 and the object therefore
      // has no level)
      if (this->present_index >= static_cast<int>(objects().n_objects()))
        this->present_index = -1;
    }
  else
    {
      while (this->present_index >=
             static_cast<int>(this->tria->levels[level()]->cells.n_objects()))
        {
          // no -> go one level up until we find
          // one with more than zero cells
          this->set_level(level() + 1);
          this->present_index = 0;
          // highest level reached?
          if (level() >= static_cast<int>(this->tria->levels.size()))
            {
              // return with past the end pointer
              this->set_level(-1);
              this->present_index = -1;
              return;
            }
        }
    }
}
 
 
template <int structdim, int dim, int spacedim>
inline void
TriaAccessorBase<structdim, dim, spacedim>::operator--()
{
  // same as operator++
  --this->present_index;
 
  if (structdim != dim)
    {
      if (this->present_index < 0)
        this->present_index = -1;
    }
  else
    {
      while (this->present_index < 0)
        {
          // no -> go one level down
          this->set_level(level() - 1);
          // lowest level reached?
          if (level() == -1)
            {
              // return with past the end pointer
              this->set_level(-1);
              this->present_index = -1;
              return;
            }
          // else
          this->present_index =
            this->tria->levels[level()]->cells.n_objects() - 1;
        }
    }
}
 
 
 
template <int structdim, int dim, int spacedim>
inline ::internal::TriangulationImplementation::TriaObjects &
TriaAccessorBase<structdim, dim, spacedim>::objects() const
{
  if (structdim == dim)
    return this->tria->levels[level()]->cells;
 
  if (structdim == 1 && dim > 1)
    return this->tria->faces->lines;
 
  if (structdim == 2 && dim > 2)
    return this->tria->faces->quads;
 
  DEAL_II_ASSERT_UNREACHABLE();
 
  return this->tria->levels[level()]->cells;
}
 
 
 
/*---------------------- Functions: InvalidAccessor -------------------------*/
 
template <int structdim, int dim, int spacedim>
InvalidAccessor<structdim, dim, spacedim>::InvalidAccessor(const void *,
                                                           const int,
                                                           const int,
                                                           const AccessorData *)
{
  Assert(false,
         ExcMessage("You are attempting an invalid conversion between "
                    "iterator/accessor types. The constructor you call "
                    "only exists to make certain template constructs "
                    "easier to write as dimension independent code but "
                    "the conversion is not valid in the current context."));
}
 
 
 
template <int structdim, int dim, int spacedim>
InvalidAccessor<structdim, dim, spacedim>::InvalidAccessor(
  const InvalidAccessor &)
{
  Assert(false,
         ExcMessage("You are attempting an invalid conversion between "
                    "iterator/accessor types. The constructor you call "
                    "only exists to make certain template constructs "
                    "easier to write as dimension independent code but "
                    "the conversion is not valid in the current context."));
}
 
 
 
template <int structdim, int dim, int spacedim>
void
InvalidAccessor<structdim, dim, spacedim>::copy_from(const InvalidAccessor &)
{
  // nothing to do here. we could throw an exception but we can't get here
  // without first creating an object which would have already thrown.
}
 
 
 
template <int structdim, int dim, int spacedim>
bool
InvalidAccessor<structdim, dim, spacedim>::operator==(
  const InvalidAccessor &) const
{
  // nothing to do here. we could throw an exception but we can't get here
  // without first creating an object which would have already thrown.
  return false;
}
 
 
 
template <int structdim, int dim, int spacedim>
bool
InvalidAccessor<structdim, dim, spacedim>::operator!=(
  const InvalidAccessor &) const
{
  // nothing to do here. we could throw an exception but we can't get here
  // without first creating an object which would have already thrown.
  return true;
}
 
 
 
template <int structdim, int dim, int spacedim>
bool
InvalidAccessor<structdim, dim, spacedim>::used() const
{
  // nothing to do here. we could throw an exception but we can't get here
  // without first creating an object which would have already thrown.
  return false;
}
 
 
 
template <int structdim, int dim, int spacedim>
bool
InvalidAccessor<structdim, dim, spacedim>::has_children() const
{
  // nothing to do here. we could throw an exception but we can't get here
  // without first creating an object which would have already thrown.
  return false;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
InvalidAccessor<structdim, dim, spacedim>::operator++() const
{}
 
 
 
template <int structdim, int dim, int spacedim>
void
InvalidAccessor<structdim, dim, spacedim>::operator--() const
{}
 
 
 
template <int structdim, int dim, int spacedim>
types::manifold_id
InvalidAccessor<structdim, dim, spacedim>::manifold_id() const
{
  return numbers::flat_manifold_id;
}
 
 
 
template <int structdim, int dim, int spacedim>
unsigned int
InvalidAccessor<structdim, dim, spacedim>::user_index() const
{
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
InvalidAccessor<structdim, dim, spacedim>::set_user_index(
  const unsigned int) const
{
  Assert(false,
         ExcMessage("You are trying to set the user index of an "
                    "invalid object."));
}
 
 
 
template <int structdim, int dim, int spacedim>
void
InvalidAccessor<structdim, dim, spacedim>::set_manifold_id(
  const types::manifold_id) const
{
  Assert(false,
         ExcMessage("You are trying to set the manifold id of an "
                    "invalid object."));
}
 
 
 
template <int structdim, int dim, int spacedim>
inline Point<spacedim> &
InvalidAccessor<structdim, dim, spacedim>::vertex(const unsigned int) const
{
  // nothing to do here. we could throw an exception but we can't get here
  // without first creating an object which would have already thrown
  static Point<spacedim> invalid_vertex;
  return invalid_vertex;
}
 
 
template <int structdim, int dim, int spacedim>
inline void *
InvalidAccessor<structdim, dim, spacedim>::line(const unsigned int) const
{
  // nothing to do here. we could throw an exception but we can't get here
  // without first creating an object which would have already thrown
  return nullptr;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline void *
InvalidAccessor<structdim, dim, spacedim>::quad(const unsigned int) const
{
  // nothing to do here. we could throw an exception but we can't get here
  // without first creating an object which would have already thrown
  return nullptr;
}
 
 
/*------------------------ Functions: TriaAccessor ---------------------------*/
 
 
namespace internal
{
  namespace TriaAccessorImplementation
  {
    // make sure that if in the following we write TriaAccessor we mean the
    // *class* ::TriaAccessor, not the enclosing namespace
    // ::internal::TriaAccessor
    using ::TriaAccessor;
 
    struct Implementation
    {
      template <int dim, int spacedim>
      static std::array<unsigned int, 1>
      get_line_indices_of_cell(const TriaAccessor<1, dim, spacedim> &)
      {
        DEAL_II_ASSERT_UNREACHABLE();
        return {};
      }
 
 
 
      template <int structdim, int dim, int spacedim>
      static std::array<unsigned int, 4>
      get_line_indices_of_cell(const TriaAccessor<2, dim, spacedim> &cell)
      {
        // For 2d cells the access cell->line_orientation() is already
        // efficient
        std::array<unsigned int, 4> line_indices = {};
        for (const unsigned int line : cell.line_indices())
          line_indices[line] = cell.line_index(line);
        return line_indices;
      }
 
      template <int structdim, int dim, int spacedim>
      static std::array<unsigned int, 12>
      get_line_indices_of_cell(
        const TriaAccessor<structdim, dim, spacedim> &cell)
      {
        std::array<unsigned int, 12> line_indices = {};
 
        // For hexahedra, the classical access via quads -> lines is too
        // inefficient. Unroll this code here to allow the compiler to inline
        // the necessary functions.
        const auto cell_reference_cell = cell.reference_cell();
        if (cell_reference_cell == ReferenceCells::Hexahedron)
          {
            constexpr auto reference_cell = ReferenceCells::Hexahedron;
            for (unsigned int f = 4; f < 6; ++f)
              {
                const auto orientation = cell.combined_face_orientation(f);
 
                // It might seem superfluous to spell out the four indices
                // that get later consumed by a for loop over these four
                // elements; however, for the compiler it is easier to inline
                // the statement of standard_to_real_face_line() when next to
                // each other, as opposed to be interleaved with a
                // line_index() call.
                const std::array<unsigned int, 4> my_indices{{
                  reference_cell.standard_to_real_face_line(0, f, orientation),
                  reference_cell.standard_to_real_face_line(1, f, orientation),
                  reference_cell.standard_to_real_face_line(2, f, orientation),
                  reference_cell.standard_to_real_face_line(3, f, orientation),
                }};
                const auto                        quad = cell.quad(f);
                for (unsigned int l = 0; l < 4; ++l)
                  line_indices[4 * (f - 4) + l] =
                    quad->line_index(my_indices[l]);
              }
            for (unsigned int f = 0; f < 2; ++f)
              {
                const auto orientation =
                  cell.get_triangulation()
                    .levels[cell.level()]
                    ->face_orientations.get_combined_orientation(cell.index(),
                                                                 f);
                const std::array<unsigned int, 2> my_indices{{
                  reference_cell.standard_to_real_face_line(0, f, orientation),
                  reference_cell.standard_to_real_face_line(1, f, orientation),
                }};
                const auto                        quad = cell.quad(f);
                line_indices[8 + f]  = quad->line_index(my_indices[0]);
                line_indices[10 + f] = quad->line_index(my_indices[1]);
              }
          }
        else if (cell_reference_cell == ReferenceCells::Tetrahedron)
          {
            constexpr auto reference_cell = ReferenceCells::Tetrahedron;
            const std::array<unsigned int, 3> orientations{
              {cell.combined_face_orientation(0),
               cell.combined_face_orientation(1),
               cell.combined_face_orientation(2)}};
            const std::array<unsigned int, 6>           my_indices{{
              reference_cell.standard_to_real_face_line(0, 0, orientations[0]),
              reference_cell.standard_to_real_face_line(1, 0, orientations[0]),
              reference_cell.standard_to_real_face_line(2, 0, orientations[0]),
              reference_cell.standard_to_real_face_line(1, 1, orientations[1]),
              reference_cell.standard_to_real_face_line(2, 1, orientations[1]),
              reference_cell.standard_to_real_face_line(1, 2, orientations[2]),
            }};
            const std::array<decltype(cell.quad(0)), 3> quads{
              {cell.quad(0), cell.quad(1), cell.quad(2)}};
            line_indices[0] = quads[0]->line_index(my_indices[0]);
            line_indices[1] = quads[0]->line_index(my_indices[1]);
            line_indices[2] = quads[0]->line_index(my_indices[2]);
            line_indices[3] = quads[1]->line_index(my_indices[3]);
            line_indices[4] = quads[1]->line_index(my_indices[4]);
            line_indices[5] = quads[2]->line_index(my_indices[5]);
          }
        else
          // For other shapes (wedges, pyramids), we do not currently
          // implement an optimized function.
          for (unsigned int l = 0; l < std::min(12U, cell.n_lines()); ++l)
            line_indices[l] = cell.line_index(l);
 
        return line_indices;
      }
 
 
 
      template <int dim, int spacedim>
      static std::pair<std::array<unsigned int, 1>,
                       std::array<types::geometric_orientation, 1>>
      get_line_indices_and_orientations_of_cell(
        const TriaAccessor<1, dim, spacedim> &)
      {
        DEAL_II_ASSERT_UNREACHABLE();
        return {};
      }
 
 
 
      template <int dim, int spacedim>
      static std::pair<std::array<unsigned int, 4>,
                       std::array<types::geometric_orientation, 4>>
      get_line_indices_and_orientations_of_cell(
        const TriaAccessor<2, dim, spacedim> &cell)
      {
        // For 2d cells the access cell->line_orientation() is already
        // efficient
        std::array<types::geometric_orientation, 4> line_orientations = {};
        std::array<unsigned int, 4>                 line_indices      = {};
        for (const unsigned int line : cell.line_indices())
          {
            line_orientations[line] = cell.line_orientation(line);
            line_indices[line]      = cell.line_index(line);
          }
        return std::make_pair(line_indices, line_orientations);
      }
 
 
 
      template <int dim, int spacedim>
      static std::pair<std::array<unsigned int, 12>,
                       std::array<types::geometric_orientation, 12>>
      get_line_indices_and_orientations_of_cell(
        const TriaAccessor<3, dim, spacedim> &cell)
      {
        std::array<types::geometric_orientation, 12> line_orientations = {};
        std::array<unsigned int, 12>                 line_indices      = {};
 
        // For hexahedra, the classical access via quads -> lines is too
        // inefficient. Unroll this code here to allow the compiler to inline
        // the necessary functions.
        const auto cell_reference_cell = cell.reference_cell();
        if (cell_reference_cell == ReferenceCells::Hexahedron)
          {
            constexpr auto reference_cell = ReferenceCells::Hexahedron;
            for (unsigned int f = 4; f < 6; ++f)
              {
                const auto orientation = cell.combined_face_orientation(f);
 
                // It might seem superfluous to spell out the four indices and
                // orientations that get later consumed by a for loop over
                // these four elements; however, for the compiler it is easier
                // to inline the statement of standard_to_real_face_line()
                // when next to each other, as opposed to be interleaved with
                // a line_index() call.
                const std::array<unsigned int, 4> my_indices{{
                  reference_cell.standard_to_real_face_line(0, f, orientation),
                  reference_cell.standard_to_real_face_line(1, f, orientation),
                  reference_cell.standard_to_real_face_line(2, f, orientation),
                  reference_cell.standard_to_real_face_line(3, f, orientation),
                }};
                const auto                        quad = cell.quad(f);
                for (unsigned int l = 0; l < 4; ++l)
                  line_indices[4 * (f - 4) + l] =
                    quad->line_index(my_indices[l]);
                const std::array<types::geometric_orientation, 4>
                  my_orientations{{
                    reference_cell.face_to_cell_line_orientation(
                      0, f, orientation, quad->line_orientation(my_indices[0])),
                    reference_cell.face_to_cell_line_orientation(
                      1, f, orientation, quad->line_orientation(my_indices[1])),
                    reference_cell.face_to_cell_line_orientation(
                      2, f, orientation, quad->line_orientation(my_indices[2])),
                    reference_cell.face_to_cell_line_orientation(
                      3, f, orientation, quad->line_orientation(my_indices[3])),
                  }};
                for (unsigned int l = 0; l < 4; ++l)
                  line_orientations[4 * (f - 4) + l] = my_orientations[l];
              }
            for (unsigned int f = 0; f < 2; ++f)
              {
                const auto orientation = cell.combined_face_orientation(f);
 
                const std::array<unsigned int, 2> my_indices{{
                  reference_cell.standard_to_real_face_line(0, f, orientation),
                  reference_cell.standard_to_real_face_line(1, f, orientation),
                }};
                const auto                        quad = cell.quad(f);
                line_indices[8 + f]  = quad->line_index(my_indices[0]);
                line_indices[10 + f] = quad->line_index(my_indices[1]);
                const std::array<types::geometric_orientation, 2>
                  my_orientations{{
                    reference_cell.face_to_cell_line_orientation(
                      0, f, orientation, quad->line_orientation(my_indices[0])),
                    reference_cell.face_to_cell_line_orientation(
                      1, f, orientation, quad->line_orientation(my_indices[1])),
                  }};
                line_orientations[8 + f]  = my_orientations[0];
                line_orientations[10 + f] = my_orientations[1];
              }
          }
        else if (cell_reference_cell == ReferenceCells::Tetrahedron)
          {
            constexpr auto reference_cell = ReferenceCells::Tetrahedron;
            const std::array<types::geometric_orientation, 3> orientations{
              {cell.combined_face_orientation(0),
               cell.combined_face_orientation(1),
               cell.combined_face_orientation(2)}};
            const std::array<unsigned int, 6>           my_indices{{
              reference_cell.standard_to_real_face_line(0, 0, orientations[0]),
              reference_cell.standard_to_real_face_line(1, 0, orientations[0]),
              reference_cell.standard_to_real_face_line(2, 0, orientations[0]),
              reference_cell.standard_to_real_face_line(1, 1, orientations[1]),
              reference_cell.standard_to_real_face_line(2, 1, orientations[1]),
              reference_cell.standard_to_real_face_line(1, 2, orientations[2]),
            }};
            const std::array<decltype(cell.quad(0)), 3> quads{
              {cell.quad(0), cell.quad(1), cell.quad(2)}};
            line_indices[0] = quads[0]->line_index(my_indices[0]);
            line_indices[1] = quads[0]->line_index(my_indices[1]);
            line_indices[2] = quads[0]->line_index(my_indices[2]);
            line_indices[3] = quads[1]->line_index(my_indices[3]);
            line_indices[4] = quads[1]->line_index(my_indices[4]);
            line_indices[5] = quads[2]->line_index(my_indices[5]);
 
            line_orientations[0] = reference_cell.face_to_cell_line_orientation(
              0, 0, orientations[0], quads[0]->line_orientation(my_indices[0]));
            line_orientations[1] = reference_cell.face_to_cell_line_orientation(
              1, 0, orientations[0], quads[0]->line_orientation(my_indices[1]));
            line_orientations[2] = reference_cell.face_to_cell_line_orientation(
              2, 0, orientations[0], quads[0]->line_orientation(my_indices[2]));
            line_orientations[3] = reference_cell.face_to_cell_line_orientation(
              1, 1, orientations[1], quads[1]->line_orientation(my_indices[3]));
            line_orientations[4] = reference_cell.face_to_cell_line_orientation(
              2, 1, orientations[1], quads[1]->line_orientation(my_indices[4]));
            line_orientations[5] = reference_cell.face_to_cell_line_orientation(
              1, 2, orientations[2], quads[2]->line_orientation(my_indices[5]));
          }
        else
          // For other shapes (wedges, pyramids), we do not currently implement
          // an optimized function
          for (unsigned int l = 0; l < std::min(12U, cell.n_lines()); ++l)
            {
              line_indices[l]      = cell.line_index(l);
              line_orientations[l] = cell.line_orientation(l);
            }
 
        return std::make_pair(line_indices, line_orientations);
      }
    };
  } // namespace TriaAccessorImplementation
} // namespace internal
 
 
 
template <int structdim, int dim, int spacedim>
inline TriaAccessor<structdim, dim, spacedim>::TriaAccessor(
  const Triangulation<dim, spacedim> *parent,
  const int                           level,
  const int                           index,
  const AccessorData                 *local_data)
  : TriaAccessorBase<structdim, dim, spacedim>(parent, level, index, local_data)
{}
 
 
 
template <int structdim, int dim, int spacedim>
inline bool
TriaAccessor<structdim, dim, spacedim>::used() const
{
  Assert(this->state() == IteratorState::valid,
         TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
           *this));
  return this->objects().used[this->present_index];
}
 
 
 
template <int structdim, int dim, int spacedim>
inline TriaIterator<TriaAccessor<0, dim, spacedim>>
TriaAccessor<structdim, dim, spacedim>::vertex_iterator(
  const unsigned int i) const
{
  return TriaIterator<TriaAccessor<0, dim, spacedim>>(this->tria,
                                                      0,
                                                      vertex_index(i));
}
 
 
 
template <int structdim, int dim, int spacedim>
inline ReferenceCell<structdim>
TriaAccessor<structdim, dim, spacedim>::reference_cell() const
{
  if constexpr (structdim == 0)
    return ReferenceCells::Vertex;
  else if constexpr (structdim == 1)
    return ReferenceCells::Line;
  else if constexpr (structdim == dim)
    return this->tria->levels[this->level()]
      ->reference_cell[this->present_index];
  else
    return this->tria->faces->get_quad_type(this->present_index);
}
 
 
 
template <int structdim, int dim, int spacedim>
inline unsigned int
TriaAccessor<structdim, dim, spacedim>::vertex_index(
  const unsigned int vertex_no) const
{
  AssertIndexRange(vertex_no, this->n_vertices());
 
  if constexpr (structdim == 1)
    {
      // This branch needs to be first (and not combined with the structdim ==
      // dim branch) so that we can get line vertex indices when setting up the
      // cell vertex index cache
      return this->objects().get_bounding_object(this->present_index,
                                                 vertex_no);
    }
  else if constexpr (structdim == dim)
    {
      // This branch should only be used after the cell vertex index cache is
      // set up
      const auto vertex_index =
        this->tria->levels[this->level()]->cached_vertex_index(
          this->present_index, vertex_no);
      Assert(
        vertex_index != numbers::invalid_unsigned_int,
        ExcMessage(
          "The present vertex's vertex_index is not valid. This is typically "
          "caused by calling cell->vertex_index() during grid refinement for "
          "a newly created cell, since the cell vertex index cache is not "
          "complete at that point."));
      return vertex_index;
    }
  else if constexpr (structdim == 2)
    {
      const auto [line_index, vertex_index] =
        this->reference_cell().standard_vertex_to_face_and_vertex_index(
          vertex_no);
      const auto vertex_within_line_index =
        this->reference_cell().standard_to_real_face_vertex(
          vertex_index, line_index, this->line_orientation(line_index));
 
      return this->line(line_index)->vertex_index(vertex_within_line_index);
    }
  else
    {
      DEAL_II_ASSERT_UNREACHABLE();
      return numbers::invalid_unsigned_int;
    }
}
 
 
 
template <int structdim, int dim, int spacedim>
inline Point<spacedim> &
TriaAccessor<structdim, dim, spacedim>::vertex(const unsigned int i) const
{
  return const_cast<Point<spacedim> &>(this->tria->vertices[vertex_index(i)]);
}
 
 
 
template <int structdim, int dim, int spacedim>
inline typename ::internal::TriangulationImplementation::
  Iterators<dim, spacedim>::line_iterator
  TriaAccessor<structdim, dim, spacedim>::line(const unsigned int i) const
{
  // checks happen in line_index
  return typename ::internal::TriangulationImplementation::
    Iterators<dim, spacedim>::line_iterator(this->tria, 0, line_index(i));
}
 
 
 
template <int structdim, int dim, int spacedim>
inline unsigned int
TriaAccessor<structdim, dim, spacedim>::line_index(const unsigned int i) const
{
  (void)i;
  AssertIndexRange(i, this->n_lines());
  Assert(structdim != 1,
         ExcMessage("You can't ask for the index of a line bounding a "
                    "one-dimensional cell because it is not bounded by "
                    "lines."));
 
  if constexpr (structdim == 2)
    {
      return this->objects().get_bounding_object(this->present_index, i);
    }
  else if constexpr (structdim == 3)
    {
      const auto [face_index, line_index] =
        this->reference_cell().standard_line_to_face_and_line_index(i);
      const auto line_within_face_index =
        this->reference_cell().standard_to_real_face_line(
          line_index, face_index, this->combined_face_orientation(face_index));
 
      return this->quad(face_index)->line_index(line_within_face_index);
    }
 
  DEAL_II_ASSERT_UNREACHABLE();
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline typename ::internal::TriangulationImplementation::
  Iterators<dim, spacedim>::quad_iterator
  TriaAccessor<structdim, dim, spacedim>::quad(const unsigned int i) const
{
  // checks happen in quad_index
  return typename ::internal::TriangulationImplementation::
    Iterators<dim, spacedim>::quad_iterator(this->tria, 0, quad_index(i));
}
 
 
 
template <int structdim, int dim, int spacedim>
inline unsigned int
TriaAccessor<structdim, dim, spacedim>::quad_index(
  const unsigned int face_no) const
{
  Assert(structdim == 3,
         ExcMessage("You can't ask for the index of a quad bounding "
                    "a one- or two-dimensional cell because it is not "
                    "bounded by quads."));
  // work around a bogus GCC-9 warning which considers face_no unused except in
  // 3d
  (void)face_no;
  if constexpr (structdim == 3)
    return this->objects().get_bounding_object(this->present_index, face_no);
  else
    return numbers::invalid_unsigned_int;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline types::geometric_orientation
TriaAccessor<structdim, dim, spacedim>::combined_face_orientation(
  const unsigned int face) const
{
  Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
  AssertIndexRange(face, n_faces());
  Assert(structdim == dim,
         ExcMessage("This function can only be used on objects "
                    "that are cells, but not on faces or edges "
                    "that bound cells."));
  // work around a bogus GCC-9 warning which considers face unused except in 3d
  (void)face;
 
  if constexpr (structdim == 1)
    return numbers::default_geometric_orientation;
  else if constexpr (structdim == 2)
    {
      // if all elements are quads (or if we have a very special consistently
      // oriented triangular mesh) then we do not store this array
      if (this->tria->levels[this->level()]->face_orientations.n_objects() == 0)
        return numbers::default_geometric_orientation;
      else
        return this->tria->levels[this->level()]
          ->face_orientations.get_combined_orientation(this->present_index,
                                                       face);
    }
  else
    return this->tria->levels[this->level()]
      ->face_orientations.get_combined_orientation(this->present_index, face);
}
 
 
 
template <int structdim, int dim, int spacedim>
inline bool
TriaAccessor<structdim, dim, spacedim>::face_orientation(
  const unsigned int face) const
{
  Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
  AssertIndexRange(face, n_faces());
  Assert(structdim == dim,
         ExcMessage("This function can only be used on objects "
                    "that are cells, but not on faces or edges "
                    "that bound cells."));
  // work around a bogus GCC-9 warning which considers face unused in 1d
  (void)face;
 
  if constexpr (structdim == 1)
    // in 1d 'faces' are vertices and those are always consistently oriented
    return true;
  else if constexpr (structdim == 2)
    return this->line_orientation(face) ==
           numbers::default_geometric_orientation;
  else
    return this->tria->levels[this->level()]->face_orientations.get_orientation(
      this->present_index, face);
}
 
 
 
template <int structdim, int dim, int spacedim>
inline bool
TriaAccessor<structdim, dim, spacedim>::face_flip(const unsigned int face) const
{
  Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
  Assert(structdim == dim,
         ExcMessage("This function can only be used on objects "
                    "that are cells, but not on faces or edges "
                    "that bound cells."));
  AssertIndexRange(face, n_faces());
  // work around a bogus GCC-9 warning which considers face unused except in 3d
  (void)face;
 
  if constexpr (structdim == 3)
    return this->tria->levels[this->level()]->face_orientations.get_flip(
      this->present_index, face);
  else
    // In 1d and 2d, face_flip is always false as faces can only be
    // 'flipped' in 3d.
    return false;
}
 
 
template <int structdim, int dim, int spacedim>
inline bool
TriaAccessor<structdim, dim, spacedim>::face_rotation(
  const unsigned int face) const
{
  Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
  Assert(structdim == dim,
         ExcMessage("This function can only be used on objects "
                    "that are cells, but not on faces or edges "
                    "that bound cells."));
  AssertIndexRange(face, n_faces());
  // work around a bogus GCC-9 warning which considers face unused except in 3d
  (void)face;
 
  if constexpr (structdim == 3)
    return this->tria->levels[this->level()]->face_orientations.get_rotation(
      this->present_index, face);
  else
    // In 1d and 2d, face_rotation is always false as faces can only be
    // 'rotated' in 3d.
    return false;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline types::geometric_orientation
TriaAccessor<structdim, dim, spacedim>::line_orientation(
  const unsigned int line) const
{
  Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
  AssertIndexRange(line, this->n_lines());
  // work around a bogus GCC-9 warning which considers line unused in 1d
  (void)line;
 
  if constexpr (structdim == 1)
    return numbers::default_geometric_orientation;
  else if constexpr (structdim == 2 && dim == 2)
    // lines in 2d are faces
    {
      const auto combined_orientation = combined_face_orientation(line);
      Assert(combined_orientation == numbers::default_geometric_orientation ||
               combined_orientation == numbers::reverse_line_orientation,
             ExcInternalError());
      return combined_orientation;
    }
  else if constexpr (structdim == 2 && dim == 3)
    {
      return this->tria->faces->get_line_orientation(this->present_index, line);
    }
  else if constexpr (structdim == 3 && dim == 3)
    {
      const auto reference_cell = this->reference_cell();
      // First pick a face on which this line is a part of, and the
      // index of the line within.
      const auto [face_index, line_index] =
        reference_cell.standard_line_to_face_and_line_index(line);
      const auto line_within_face_index =
        reference_cell.standard_to_real_face_line(
          line_index, face_index, this->combined_face_orientation(face_index));
 
      // Then query how that line is oriented within that face:
      return reference_cell.face_to_cell_line_orientation(
        line_index,
        face_index,
        this->combined_face_orientation(face_index),
        this->quad(face_index)->line_orientation(line_within_face_index));
    }
  else
    {
      DEAL_II_ASSERT_UNREACHABLE();
      return numbers::invalid_geometric_orientation;
    }
}
 
 
 
template <int structdim, int dim, int spacedim>
inline void
TriaAccessor<structdim, dim, spacedim>::set_line_orientation(
  const unsigned int                 line,
  const types::geometric_orientation value) const
{
  Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
  AssertIndexRange(line, this->n_lines());
  Assert(dim != 1,
         ExcMessage("In 1d lines are cells and thus do not need to have their "
                    "orientations set."));
  Assert(dim != 2,
         ExcMessage("In 2d lines are faces, and, for compatibility with other "
                    "dimensions, their orientations should be set via "
                    "set_combined_face_orientation()."));
  // work around a bogus GCC-9 warning which considers line and value unused
  // except in 3d
  (void)line;
  (void)value;
 
  if constexpr (dim == 3)
    {
      // We set line orientations per face, not per cell, so this only works for
      // faces in 3d.
      Assert(structdim == 2, ExcNotImplemented());
      this->tria->faces->set_line_orientation(this->present_index, line, value);
    }
}
 
 
 
template <int structdim, int dim, int spacedim>
inline void
TriaAccessor<structdim, dim, spacedim>::set_combined_face_orientation(
  const unsigned int                 face_no,
  const types::geometric_orientation combined_orientation) const
{
  Assert(used(), TriaAccessorExceptions::ExcCellNotUsed());
  AssertIndexRange(face_no, this->n_faces());
  Assert(structdim == dim,
         ExcMessage("This function can only be used on objects that are cells "
                    "and not on objects which bound cells."));
  AssertIndexRange(combined_orientation,
                   reference_cell().n_face_orientations(face_no));
 
  // face_orientations is not set up in 1d
  if (dim > 1)
    this->tria->levels[this->level()]
      ->face_orientations.set_combined_orientation(this->present_index,
                                                   face_no,
                                                   combined_orientation);
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::set_used_flag() const
{
  Assert(this->state() == IteratorState::valid,
         TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
           *this));
  this->objects().used[this->present_index] = true;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::clear_used_flag() const
{
  Assert(this->state() == IteratorState::valid,
         TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
           *this));
  this->objects().used[this->present_index] = false;
}
 
 
template <int structdim, int dim, int spacedim>
int
TriaAccessor<structdim, dim, spacedim>::child_index(const unsigned int i) const
{
  Assert(has_children(), TriaAccessorExceptions::ExcCellHasNoChildren());
  AssertIndexRange(i, n_children());
 
  // each set of two children are stored
  // consecutively, so we only have to find
  // the location of the set of children
  const auto        &objects       = this->objects();
  const unsigned int n_sets_of_two = objects.children_per_object / 2;
  return objects.children[n_sets_of_two * this->present_index + i / 2] + i % 2;
}
 
 
 
template <int structdim, int dim, int spacedim>
int
TriaAccessor<structdim, dim, spacedim>::isotropic_child_index(
  const unsigned int i) const
{
  AssertIndexRange(i, this->objects().children_per_object);
 
  switch (structdim)
    {
      case 1:
        return child_index(i);
      case 2:
        {
          const RefinementCase<2> this_refinement_case(
            static_cast<std::uint8_t>(refinement_case()));
 
          Assert(this_refinement_case != RefinementCase<2>::no_refinement,
                 TriaAccessorExceptions::ExcCellHasNoChildren());
 
          if (this_refinement_case == RefinementCase<2>::cut_xy)
            return child_index(i);
          else if ((this_refinement_case == RefinementCase<2>::cut_x) &&
                   (child(i % 2)->refinement_case() ==
                    RefinementCase<2>::cut_y))
            return child(i % 2)->child_index(i / 2);
          else if ((this_refinement_case == RefinementCase<2>::cut_y) &&
                   (child(i / 2)->refinement_case() ==
                    RefinementCase<2>::cut_x))
            return child(i / 2)->child_index(i % 2);
          else
            Assert(
              false,
              ExcMessage(
                "This cell has no grandchildren equivalent to isotropic refinement"));
          break;
        }
 
      case 3:
        DEAL_II_NOT_IMPLEMENTED();
    }
  return -1;
}
 
 
 
template <int structdim, int dim, int spacedim>
RefinementCase<structdim>
TriaAccessor<structdim, dim, spacedim>::refinement_case() const
{
  Assert(this->state() == IteratorState::valid,
         TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
           *this));
 
  switch (structdim)
    {
      case 1:
        return (RefinementCase<structdim>(
          this->objects().children[this->present_index] != -1 ?
            // cast the branches here first to uchar and then (above) to
            // RefinementCase<structdim> so that the conversion is valid even
            // for the case structdim>1 (for which this part of the code is dead
            // anyway)
            static_cast<std::uint8_t>(RefinementCase<1>::cut_x) :
            static_cast<std::uint8_t>(RefinementCase<1>::no_refinement)));
 
      default:
        Assert(static_cast<unsigned int>(this->present_index) <
                 this->objects().refinement_cases.size(),
               ExcIndexRange(this->present_index,
                             0,
                             this->objects().refinement_cases.size()));
 
        return (static_cast<RefinementCase<structdim>>(
          this->objects().refinement_cases[this->present_index]));
    }
}
 
 
 
template <int structdim, int dim, int spacedim>
inline TriaIterator<TriaAccessor<structdim, dim, spacedim>>
TriaAccessor<structdim, dim, spacedim>::child(const unsigned int i) const
 
{
  // checking of 'i' happens in child_index
  const TriaIterator<TriaAccessor<structdim, dim, spacedim>> q(
    this->tria, (dim == structdim ? this->level() + 1 : 0), child_index(i));
 
  Assert((q.state() == IteratorState::past_the_end) || q->used(),
         ExcInternalError());
 
  return q;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline unsigned int
TriaAccessor<structdim, dim, spacedim>::child_iterator_to_index(
  const TriaIterator<TriaAccessor<structdim, dim, spacedim>> &child) const
{
  const auto n_children = this->n_children();
  for (unsigned int child_n = 0; child_n < n_children; ++child_n)
    if (this->child(child_n) == child)
      return child_n;
 
  Assert(false,
         ExcMessage("The given child is not a child of the current object."));
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline TriaIterator<TriaAccessor<structdim, dim, spacedim>>
TriaAccessor<structdim, dim, spacedim>::isotropic_child(
  const unsigned int i) const
{
  // checking of 'i' happens in child() or child_index() called below
  switch (structdim)
    {
      case 1:
        // no anisotropic refinement in 1d
        return child(i);
 
      case 2:
        {
          const RefinementCase<2> this_refinement_case(
            static_cast<std::uint8_t>(refinement_case()));
 
          Assert(this_refinement_case != RefinementCase<2>::no_refinement,
                 TriaAccessorExceptions::ExcCellHasNoChildren());
 
          if (this_refinement_case == RefinementCase<2>::cut_xy)
            return child(i);
          else if ((this_refinement_case == RefinementCase<2>::cut_x) &&
                   (child(i % 2)->refinement_case() ==
                    RefinementCase<2>::cut_y))
            return child(i % 2)->child(i / 2);
          else if ((this_refinement_case == RefinementCase<2>::cut_y) &&
                   (child(i / 2)->refinement_case() ==
                    RefinementCase<2>::cut_x))
            return child(i / 2)->child(i % 2);
          else
            Assert(
              false,
              ExcMessage(
                "This cell has no grandchildren equivalent to isotropic refinement"));
          break;
        }
 
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
  // we don't get here but have to return
  // something...
  return child(0);
}
 
 
 
template <int structdim, int dim, int spacedim>
inline bool
TriaAccessor<structdim, dim, spacedim>::has_children() const
{
  Assert(this->state() == IteratorState::valid,
         TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
           *this));
 
  // each set of two children are stored consecutively, so we only have to find
  // the location of the set of children
  const auto        &objects       = this->objects();
  const unsigned int n_sets_of_two = objects.children_per_object / 2;
  return (objects.children[n_sets_of_two * this->present_index] != -1);
}
 
 
 
template <int structdim, int dim, int spacedim>
inline unsigned int
TriaAccessor<structdim, dim, spacedim>::n_children() const
{
  Assert(this->state() == IteratorState::valid,
         TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
           *this));
 
  return reference_cell().n_children(refinement_case());
}
 
 
 
template <int structdim, int dim, int spacedim>
inline void
TriaAccessor<structdim, dim, spacedim>::set_refinement_case(
  const RefinementCase<structdim> &refinement_case) const
{
  Assert(this->state() == IteratorState::valid,
         TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
           *this));
  Assert(static_cast<unsigned int>(this->present_index) <
           this->objects().refinement_cases.size(),
         ExcIndexRange(this->present_index,
                       0,
                       this->objects().refinement_cases.size()));
 
  this->objects().refinement_cases[this->present_index] = refinement_case;
}
 
 
template <int structdim, int dim, int spacedim>
inline void
TriaAccessor<structdim, dim, spacedim>::clear_refinement_case() const
{
  Assert(this->state() == IteratorState::valid,
         TriaAccessorExceptions::ExcDereferenceInvalidObject<TriaAccessor>(
           *this));
  Assert(static_cast<unsigned int>(this->present_index) <
           this->objects().refinement_cases.size(),
         ExcIndexRange(this->present_index,
                       0,
                       this->objects().refinement_cases.size()));
 
  this->objects().refinement_cases[this->present_index] =
    RefinementCase<structdim>::no_refinement;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::set_children(const unsigned int i,
                                                     const int index) const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  Assert(i % 2 == 0, TriaAccessorExceptions::ExcSetOnlyEvenChildren(i));
 
  // each set of two children are stored consecutively, so we only have to find
  // the location of the set of children
  auto              &objects       = this->objects();
  const unsigned int n_sets_of_two = objects.children_per_object / 2;
  Assert(
    // clearing the child index for a cell
    (index == -1) ||
      // if setting the child index for the i'th child (with i==0), then the
      // index must be a non-negative number
      (i == 0 && !this->has_children() && (index >= 0)) ||
      // if setting the child index for the i'th child (with i>0), then the
      // previously stored index must be the invalid index
      (i > 0 && this->has_children() && (index >= 0) &&
       objects.children[n_sets_of_two * this->present_index + i / 2] == -1),
    TriaAccessorExceptions::ExcCantSetChildren(index));
 
  objects.children[n_sets_of_two * this->present_index + i / 2] = index;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::clear_children() const
{
  // each set of two children are stored consecutively, so we only have to find
  // the location of the set of children
  const auto        &objects       = this->objects();
  const unsigned int n_sets_of_two = objects.children_per_object / 2;
  for (unsigned int i = 0; i < n_sets_of_two; ++i)
    set_children(2 * i, -1);
}
 
 
 
template <int structdim, int dim, int spacedim>
inline bool
TriaAccessor<structdim, dim, spacedim>::user_flag_set() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  return this->objects().user_flags[this->present_index];
}
 
 
 
template <int structdim, int dim, int spacedim>
inline void
TriaAccessor<structdim, dim, spacedim>::set_user_flag() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  this->objects().user_flags[this->present_index] = true;
}
 
 
 
template <int structdim, int dim, int spacedim>
inline void
TriaAccessor<structdim, dim, spacedim>::clear_user_flag() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  this->objects().user_flags[this->present_index] = false;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::recursively_set_user_flag() const
{
  set_user_flag();
 
  if (this->has_children())
    for (unsigned int c = 0; c < this->n_children(); ++c)
      this->child(c)->recursively_set_user_flag();
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::recursively_clear_user_flag() const
{
  clear_user_flag();
 
  if (this->has_children())
    for (unsigned int c = 0; c < this->n_children(); ++c)
      this->child(c)->recursively_clear_user_flag();
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::clear_user_data() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  this->objects().clear_user_data(this->present_index);
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::set_user_pointer(void *p) const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  this->objects().user_pointer(this->present_index) = p;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::clear_user_pointer() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  this->objects().user_pointer(this->present_index) = nullptr;
}
 
 
 
template <int structdim, int dim, int spacedim>
void *
TriaAccessor<structdim, dim, spacedim>::user_pointer() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  return this->objects().user_pointer(this->present_index);
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::recursively_set_user_pointer(
  void *p) const
{
  set_user_pointer(p);
 
  if (this->has_children())
    for (unsigned int c = 0; c < this->n_children(); ++c)
      this->child(c)->recursively_set_user_pointer(p);
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::recursively_clear_user_pointer() const
{
  clear_user_pointer();
 
  if (this->has_children())
    for (unsigned int c = 0; c < this->n_children(); ++c)
      this->child(c)->recursively_clear_user_pointer();
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::set_user_index(
  const unsigned int p) const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  this->objects().user_index(this->present_index) = p;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::clear_user_index() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  this->objects().user_index(this->present_index) = 0;
}
 
 
 
template <int structdim, int dim, int spacedim>
unsigned int
TriaAccessor<structdim, dim, spacedim>::user_index() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  return this->objects().user_index(this->present_index);
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::recursively_set_user_index(
  const unsigned int p) const
{
  set_user_index(p);
 
  if (this->has_children())
    for (unsigned int c = 0; c < this->n_children(); ++c)
      this->child(c)->recursively_set_user_index(p);
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::recursively_clear_user_index() const
{
  clear_user_index();
 
  if (this->has_children())
    for (unsigned int c = 0; c < this->n_children(); ++c)
      this->child(c)->recursively_clear_user_index();
}
 
 
 
template <int structdim, int dim, int spacedim>
inline unsigned int
TriaAccessor<structdim, dim, spacedim>::max_refinement_depth() const
{
  if (!this->has_children())
    return 0;
 
  unsigned int max_depth = 1;
  for (unsigned int c = 0; c < n_children(); ++c)
    max_depth = std::max(max_depth, child(c)->max_refinement_depth() + 1);
  return max_depth;
}
 
 
 
template <int structdim, int dim, int spacedim>
unsigned int
TriaAccessor<structdim, dim, spacedim>::n_active_descendants() const
{
  if (!this->has_children())
    return 1;
  else
    {
      unsigned int sum = 0;
      for (unsigned int c = 0; c < n_children(); ++c)
        sum += this->child(c)->n_active_descendants();
      return sum;
    }
}
 
 
 
template <int structdim, int dim, int spacedim>
types::boundary_id
TriaAccessor<structdim, dim, spacedim>::boundary_id() const
{
  Assert(structdim < dim, ExcImpossibleInDim(dim));
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
 
  return this->objects()
    .boundary_or_material_id[this->present_index]
    .boundary_id;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::set_boundary_id(
  const types::boundary_id boundary_ind) const
{
  Assert(structdim < dim, ExcImpossibleInDim(dim));
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  Assert(boundary_ind != numbers::internal_face_boundary_id,
         ExcMessage("You are trying to set the boundary_id to an invalid "
                    "value (numbers::internal_face_boundary_id is reserved)."));
  Assert(this->at_boundary(),
         ExcMessage("You are trying to set the boundary_id of an "
                    "internal object, which is not allowed!"));
 
  this->objects().boundary_or_material_id[this->present_index].boundary_id =
    boundary_ind;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::set_boundary_id_internal(
  const types::boundary_id boundary_ind) const
{
  Assert(structdim < dim, ExcImpossibleInDim(dim));
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
 
  this->objects().boundary_or_material_id[this->present_index].boundary_id =
    boundary_ind;
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::set_all_boundary_ids(
  const types::boundary_id boundary_ind) const
{
  set_boundary_id(boundary_ind);
 
  switch (structdim)
    {
      case 1:
        // 1d objects have no sub-objects where we have to do anything
        break;
 
      case 2:
        // for boundary quads also set boundary_id of bounding lines
        for (unsigned int i = 0; i < this->n_lines(); ++i)
          this->line(i)->set_boundary_id(boundary_ind);
        break;
 
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
}
 
 
 
template <int structdim, int dim, int spacedim>
bool
TriaAccessor<structdim, dim, spacedim>::at_boundary() const
{
  // error checking is done in boundary_id()
  return (boundary_id() != numbers::internal_face_boundary_id);
}
 
 
 
template <int structdim, int dim, int spacedim>
const Manifold<dim, spacedim> &
TriaAccessor<structdim, dim, spacedim>::get_manifold() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  return this->tria->get_manifold(this->manifold_id());
}
 
 
template <int structdim, int dim, int spacedim>
types::manifold_id
TriaAccessor<structdim, dim, spacedim>::manifold_id() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
 
  return this->objects().manifold_id[this->present_index];
}
 
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::set_manifold_id(
  const types::manifold_id manifold_ind) const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
 
  this->objects().manifold_id[this->present_index] = manifold_ind;
}
 
 
template <int structdim, int dim, int spacedim>
void
TriaAccessor<structdim, dim, spacedim>::set_all_manifold_ids(
  const types::manifold_id manifold_ind) const
{
  set_manifold_id(manifold_ind);
 
  if (this->has_children())
    for (unsigned int c = 0; c < this->n_children(); ++c)
      this->child(c)->set_all_manifold_ids(manifold_ind);
 
  switch (structdim)
    {
      case 1:
        if (dim == 1)
          {
            (*this->tria->vertex_to_manifold_id_map_1d)[vertex_index(0)] =
              manifold_ind;
            (*this->tria->vertex_to_manifold_id_map_1d)[vertex_index(1)] =
              manifold_ind;
          }
        break;
 
      case 2:
        // for quads/simplices also set manifold_id of bounding lines
        for (unsigned int i = 0; i < this->n_lines(); ++i)
          this->line(i)->set_manifold_id(manifold_ind);
        break;
      default:
        DEAL_II_NOT_IMPLEMENTED();
    }
}
 
 
 
template <int structdim, int dim, int spacedim>
double
TriaAccessor<structdim, dim, spacedim>::diameter() const
{
  std_cxx26::inplace_vector<Point<spacedim>,
#  ifndef _MSC_VER
                            ReferenceCells::max_n_vertices<structdim>()
#  else
                            GeometryInfo<structdim>::vertices_per_cell
#  endif
                            >
    vertices(this->n_vertices());
 
  for (unsigned int v = 0; v < vertices.size(); ++v)
    vertices[v] = this->vertex(v);
 
  return internal::TriaAccessorImplementation::diameter<structdim, spacedim>(
    vertices);
}
 
 
 
template <int dim, int spacedim>
double
CellAccessor<dim, spacedim>::diameter(
  const Mapping<dim, spacedim> &mapping) const
{
  return internal::TriaAccessorImplementation::diameter<dim, spacedim>(
    mapping.get_vertices(typename Triangulation<dim, spacedim>::cell_iterator(
      this->tria, this->level(), this->index())));
}
 
 
 
template <int structdim, int dim, int spacedim>
std::pair<Point<spacedim>, double>
TriaAccessor<structdim, dim, spacedim>::enclosing_ball() const
{
  // If the object is one dimensional, the enclosing ball is the initial iterate
  // i.e., the ball's center and diameter are the center and the diameter of the
  // object.
  if (structdim == 1)
    return std::make_pair((this->vertex(1) + this->vertex(0)) * 0.5,
                          (this->vertex(1) - this->vertex(0)).norm() * 0.5);
 
  // The list is_initial_guess_vertex contains bool values and has the same size
  // as the number of vertices per object. The entries of
  // is_initial_guess_vertex are set true only for those two vertices
  // corresponding to the largest diagonal which is being used to construct the
  // initial ball. We employ this mask to skip these two vertices while
  // enlarging the ball.
  std::vector<bool> is_initial_guess_vertex(this->n_vertices());
 
  // First let all the vertices be outside
  std::fill(is_initial_guess_vertex.begin(),
            is_initial_guess_vertex.end(),
            false);
 
  // Get an initial guess by looking at the largest diagonal
  Point<spacedim> center;
  double          radius = 0;
 
  switch (structdim)
    {
      case 2:
        {
          const Point<spacedim> p30(this->vertex(3) - this->vertex(0));
          const Point<spacedim> p21(this->vertex(2) - this->vertex(1));
          if (p30.norm() > p21.norm())
            {
              center                     = this->vertex(0) + 0.5 * p30;
              radius                     = p30.norm() / 2.;
              is_initial_guess_vertex[3] = true;
              is_initial_guess_vertex[0] = true;
            }
          else
            {
              center                     = this->vertex(1) + 0.5 * p21;
              radius                     = p21.norm() / 2.;
              is_initial_guess_vertex[2] = true;
              is_initial_guess_vertex[1] = true;
            }
          break;
        }
      case 3:
        {
          const Point<spacedim>     p70(this->vertex(7) - this->vertex(0));
          const Point<spacedim>     p61(this->vertex(6) - this->vertex(1));
          const Point<spacedim>     p25(this->vertex(2) - this->vertex(5));
          const Point<spacedim>     p34(this->vertex(3) - this->vertex(4));
          const std::vector<double> diagonals = {p70.norm(),
                                                 p61.norm(),
                                                 p25.norm(),
                                                 p34.norm()};
          const std::vector<double>::const_iterator it =
            std::max_element(diagonals.begin(), diagonals.end());
          if (it == diagonals.begin())
            {
              center                     = this->vertex(0) + 0.5 * p70;
              is_initial_guess_vertex[7] = true;
              is_initial_guess_vertex[0] = true;
            }
          else if (it == diagonals.begin() + 1)
            {
              center                     = this->vertex(1) + 0.5 * p61;
              is_initial_guess_vertex[6] = true;
              is_initial_guess_vertex[1] = true;
            }
          else if (it == diagonals.begin() + 2)
            {
              center                     = this->vertex(5) + 0.5 * p25;
              is_initial_guess_vertex[2] = true;
              is_initial_guess_vertex[5] = true;
            }
          else
            {
              center                     = this->vertex(4) + 0.5 * p34;
              is_initial_guess_vertex[3] = true;
              is_initial_guess_vertex[4] = true;
            }
          radius = *it * 0.5;
          break;
        }
      default:
        DEAL_II_NOT_IMPLEMENTED();
        return std::pair<Point<spacedim>, double>();
    }
 
  // For each vertex that is found to be geometrically outside the ball
  // enlarge the ball  so that the new ball contains both the previous ball
  // and the given vertex.
  for (const unsigned int v : this->vertex_indices())
    if (!is_initial_guess_vertex[v])
      {
        const double distance = center.distance(this->vertex(v));
        if (distance > radius)
          {
            // we found a vertex which is outside of the ball extend it (move
            // center and change radius)
            const Point<spacedim> pCV(center - this->vertex(v));
            radius = (distance + radius) * 0.5;
            center = this->vertex(v) + pCV * (radius / distance);
 
            // Now the new ball constructed in this block encloses the vertex
            // (v) that was found to be geometrically outside the old ball.
          }
      }
  if constexpr (running_in_debug_mode())
    {
      bool all_vertices_within_ball = true;
 
      // Set all_vertices_within_ball false if any of the vertices of the object
      // are geometrically outside the ball
      for (const unsigned int v : this->vertex_indices())
        if (center.distance(this->vertex(v)) >
            radius + 100. * std::numeric_limits<double>::epsilon())
          {
            all_vertices_within_ball = false;
            break;
          }
      // If all the vertices are not within the ball throw error
      Assert(all_vertices_within_ball, ExcInternalError());
    }
  return std::make_pair(center, radius);
}
 
 
template <int structdim, int dim, int spacedim>
double
TriaAccessor<structdim, dim, spacedim>::minimum_vertex_distance() const
{
  switch (structdim)
    {
      case 1:
        return (this->vertex(1) - this->vertex(0)).norm();
      case 2:
      case 3:
        {
          double min = std::numeric_limits<double>::max();
          for (const unsigned int i : this->vertex_indices())
            for (unsigned int j = i + 1; j < this->n_vertices(); ++j)
              min = std::min(min,
                             (this->vertex(i) - this->vertex(j)) *
                               (this->vertex(i) - this->vertex(j)));
          return std::sqrt(min);
        }
      default:
        DEAL_II_NOT_IMPLEMENTED();
        return -1e10;
    }
}
 
 
template <int structdim, int dim, int spacedim>
bool
TriaAccessor<structdim, dim, spacedim>::is_translation_of(
  const TriaIterator<TriaAccessor<structdim, dim, spacedim>> &o) const
{
  // go through the vertices and check... The cell is a translation of the
  // previous one in case the distance between the individual vertices in the
  // two cell is the same for all the vertices. So do the check by first getting
  // the distance on the first vertex, and then checking whether all others have
  // the same down to rounding errors (we have to be careful here because the
  // calculation of the displacement between one cell and the next can already
  // result in the loss of one or two digits), so we choose 1e-12 times the
  // distance between the zeroth vertices here.
  bool                      is_translation = true;
  const Tensor<1, spacedim> dist           = o->vertex(0) - this->vertex(0);
  const double              tol_square     = 1e-24 * dist.norm_square();
  for (unsigned int i = 1; i < this->n_vertices(); ++i)
    {
      const Tensor<1, spacedim> dist_new =
        (o->vertex(i) - this->vertex(i)) - dist;
      if (dist_new.norm_square() > tol_square)
        {
          is_translation = false;
          break;
        }
    }
  return is_translation;
}
 
 
 
template <int structdim, int dim, int spacedim>
unsigned int
TriaAccessor<structdim, dim, spacedim>::n_vertices() const
{
  return this->reference_cell().n_vertices();
}
 
 
 
template <int structdim, int dim, int spacedim>
unsigned int
TriaAccessor<structdim, dim, spacedim>::n_lines() const
{
  return this->reference_cell().n_lines();
}
 
 
 
template <int structdim, int dim, int spacedim>
unsigned int
TriaAccessor<structdim, dim, spacedim>::n_faces() const
{
  Assert(structdim == dim,
         ExcMessage("This function can only be used on objects "
                    "that are cells, but not on faces or edges "
                    "that bound cells."));
 
  return this->reference_cell().n_faces();
}
 
 
 
template <int structdim, int dim, int spacedim>
std_cxx20::ranges::iota_view<unsigned int, unsigned int>
TriaAccessor<structdim, dim, spacedim>::vertex_indices() const
{
  return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
                                                                  n_vertices());
}
 
 
 
template <int structdim, int dim, int spacedim>
std_cxx20::ranges::iota_view<unsigned int, unsigned int>
TriaAccessor<structdim, dim, spacedim>::line_indices() const
{
  return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
                                                                  n_lines());
}
 
 
 
template <int structdim, int dim, int spacedim>
std_cxx20::ranges::iota_view<unsigned int, unsigned int>
TriaAccessor<structdim, dim, spacedim>::face_indices() const
{
  return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
                                                                  n_faces());
}
 
 
 
/*----------------- Functions: TriaAccessor<0,dim,spacedim> -----------------*/
 
template <int dim, int spacedim>
inline TriaAccessor<0, dim, spacedim>::TriaAccessor(
  const Triangulation<dim, spacedim> *tria,
  const unsigned int                  vertex_index)
  : tria(tria)
  , global_vertex_index(vertex_index)
{}
 
 
 
template <int dim, int spacedim>
inline TriaAccessor<0, dim, spacedim>::TriaAccessor(
  const Triangulation<dim, spacedim> *tria,
  const int /*level*/,
  const int index,
  const AccessorData *)
  : tria(tria)
  , global_vertex_index(index)
{}
 
 
 
template <int dim, int spacedim>
template <int structdim2, int dim2, int spacedim2>
inline TriaAccessor<0, dim, spacedim>::TriaAccessor(
  const TriaAccessor<structdim2, dim2, spacedim2> &)
  : tria(nullptr)
  , global_vertex_index(numbers::invalid_unsigned_int)
{
  Assert(false, ExcImpossibleInDim(0));
}
 
 
 
template <int dim, int spacedim>
template <int structdim2, int dim2, int spacedim2>
inline TriaAccessor<0, dim, spacedim>::TriaAccessor(
  const InvalidAccessor<structdim2, dim2, spacedim2> &)
  : tria(nullptr)
  , global_vertex_index(numbers::invalid_unsigned_int)
{
  Assert(false, ExcImpossibleInDim(0));
}
 
 
 
template <int dim, int spacedim>
inline void
TriaAccessor<0, dim, spacedim>::copy_from(const TriaAccessor &t)
{
  tria                = t.tria;
  global_vertex_index = t.global_vertex_index;
}
 
 
 
template <int dim, int spacedim>
inline bool
TriaAccessor<0, dim, spacedim>::operator<(
  const TriaAccessor<0, dim, spacedim> &other) const
{
  Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators());
 
  return (global_vertex_index < other.global_vertex_index);
}
 
 
 
template <int dim, int spacedim>
inline IteratorState::IteratorStates
TriaAccessor<0, dim, spacedim>::state() const
{
  if (global_vertex_index != numbers::invalid_unsigned_int)
    return IteratorState::valid;
  else
    return IteratorState::past_the_end;
}
 
 
 
template <int dim, int spacedim>
inline int
TriaAccessor<0, dim, spacedim>::level()
{
  return 0;
}
 
 
 
template <int dim, int spacedim>
inline int
TriaAccessor<0, dim, spacedim>::index() const
{
  return global_vertex_index;
}
 
 
 
template <int dim, int spacedim>
inline const Triangulation<dim, spacedim> &
TriaAccessor<0, dim, spacedim>::get_triangulation() const
{
  return *tria;
}
 
 
 
template <int dim, int spacedim>
inline void
TriaAccessor<0, dim, spacedim>::operator++()
{
  ++global_vertex_index;
  if (global_vertex_index >= tria->n_vertices())
    global_vertex_index = numbers::invalid_unsigned_int;
}
 
 
 
template <int dim, int spacedim>
inline void
TriaAccessor<0, dim, spacedim>::operator--()
{
  if (global_vertex_index != numbers::invalid_unsigned_int)
    {
      if (global_vertex_index != 0)
        --global_vertex_index;
      else
        global_vertex_index = numbers::invalid_unsigned_int;
    }
}
 
 
 
template <int dim, int spacedim>
inline bool
TriaAccessor<0, dim, spacedim>::operator==(const TriaAccessor &t) const
{
  const bool result =
    ((tria == t.tria) && (global_vertex_index == t.global_vertex_index));
 
  return result;
}
 
 
 
template <int dim, int spacedim>
inline bool
TriaAccessor<0, dim, spacedim>::operator!=(const TriaAccessor &t) const
{
  return !(*this == t);
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
TriaAccessor<0, dim, spacedim>::vertex_index(const unsigned int) const
{
  return global_vertex_index;
}
 
 
 
template <int dim, int spacedim>
inline Point<spacedim> &
TriaAccessor<0, dim, spacedim>::vertex(const unsigned int) const
{
  return const_cast<Point<spacedim> &>(
    this->tria->vertices[global_vertex_index]);
}
 
 
 
template <int dim, int spacedim>
inline typename ::internal::TriangulationImplementation::
  Iterators<dim, spacedim>::line_iterator
  TriaAccessor<0, dim, spacedim>::line(const unsigned int)
{
  return typename ::internal::TriangulationImplementation::
    Iterators<dim, spacedim>::line_iterator();
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
TriaAccessor<0, dim, spacedim>::line_index(const unsigned int)
{
  Assert(false, ExcImpossibleInDim(0));
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim, int spacedim>
inline typename ::internal::TriangulationImplementation::
  Iterators<dim, spacedim>::quad_iterator
  TriaAccessor<0, dim, spacedim>::quad(const unsigned int)
{
  return typename ::internal::TriangulationImplementation::
    Iterators<dim, spacedim>::quad_iterator();
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
TriaAccessor<0, dim, spacedim>::quad_index(const unsigned int)
{
  Assert(false, ExcImpossibleInDim(0));
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim, int spacedim>
inline double
TriaAccessor<0, dim, spacedim>::diameter() const
{
  return 0.;
}
 
 
 
template <int dim, int spacedim>
inline double
TriaAccessor<0, dim, spacedim>::extent_in_direction(const unsigned int) const
{
  return 0.;
}
 
 
 
template <int dim, int spacedim>
inline Point<spacedim>
TriaAccessor<0, dim, spacedim>::center(const bool, const bool) const
{
  return this->tria->vertices[global_vertex_index];
}
 
 
 
template <int dim, int spacedim>
inline double
TriaAccessor<0, dim, spacedim>::measure() const
{
  return 1.0;
}
 
 
 
template <int dim, int spacedim>
inline types::geometric_orientation
TriaAccessor<0, dim, spacedim>::combined_face_orientation(
  const unsigned int /*face*/)
{
  return numbers::default_geometric_orientation;
}
 
 
 
template <int dim, int spacedim>
inline bool
TriaAccessor<0, dim, spacedim>::face_orientation(const unsigned int /*face*/)
{
  return false;
}
 
 
 
template <int dim, int spacedim>
inline bool
TriaAccessor<0, dim, spacedim>::face_flip(const unsigned int /*face*/)
{
  return false;
}
 
 
 
template <int dim, int spacedim>
inline bool
TriaAccessor<0, dim, spacedim>::face_rotation(const unsigned int /*face*/)
{
  return false;
}
 
 
 
template <int dim, int spacedim>
inline types::geometric_orientation
TriaAccessor<0, dim, spacedim>::line_orientation(const unsigned int /*line*/)
{
  return numbers::reverse_line_orientation;
}
 
 
 
template <int dim, int spacedim>
inline bool
TriaAccessor<0, dim, spacedim>::has_children()
{
  return false;
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
TriaAccessor<0, dim, spacedim>::n_children()
{
  return 0;
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
TriaAccessor<0, dim, spacedim>::n_active_descendants()
{
  return 0;
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
TriaAccessor<0, dim, spacedim>::max_refinement_depth()
{
  return 0;
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
TriaAccessor<0, dim, spacedim>::child_iterator_to_index(
  const TriaIterator<TriaAccessor<0, dim, spacedim>> &)
{
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim, int spacedim>
inline TriaIterator<TriaAccessor<0, dim, spacedim>>
TriaAccessor<0, dim, spacedim>::child(const unsigned int)
{
  return TriaIterator<TriaAccessor<0, dim, spacedim>>();
}
 
 
 
template <int dim, int spacedim>
inline TriaIterator<TriaAccessor<0, dim, spacedim>>
TriaAccessor<0, dim, spacedim>::isotropic_child(const unsigned int)
{
  return TriaIterator<TriaAccessor<0, dim, spacedim>>();
}
 
 
 
template <int dim, int spacedim>
inline RefinementCase<0>
TriaAccessor<0, dim, spacedim>::refinement_case()
{
  return {RefinementPossibilities<0>::no_refinement};
}
 
 
 
template <int dim, int spacedim>
inline int
TriaAccessor<0, dim, spacedim>::child_index(const unsigned int)
{
  return -1;
}
 
 
 
template <int dim, int spacedim>
inline int
TriaAccessor<0, dim, spacedim>::isotropic_child_index(const unsigned int)
{
  return -1;
}
 
 
 
template <int dim, int spacedim>
inline bool
TriaAccessor<0, dim, spacedim>::used() const
{
  return tria->vertex_used(global_vertex_index);
}
 
 
 
/*------------------- Functions: TriaAccessor<0,1,spacedim> -----------------*/
 
template <int spacedim>
inline TriaAccessor<0, 1, spacedim>::TriaAccessor(
  const Triangulation<1, spacedim> *tria,
  const VertexKind                  vertex_kind,
  const unsigned int                vertex_index)
  : tria(tria)
  , vertex_kind(vertex_kind)
  , global_vertex_index(vertex_index)
{}
 
 
 
template <int spacedim>
inline TriaAccessor<0, 1, spacedim>::TriaAccessor(
  const Triangulation<1, spacedim> *tria,
  const int                         level,
  const int                         index,
  const AccessorData *)
  : tria(tria)
  , vertex_kind(interior_vertex)
  , global_vertex_index(numbers::invalid_unsigned_int)
{
  // in general, calling this constructor should yield an error -- users should
  // instead call the one immediately above. however, if you create something
  // like Triangulation<1>::face_iterator() then this calls the default
  // constructor of the iterator which calls the accessor with argument list
  // (0,-2,-2,0), so in this particular case accept this call and create an
  // object that corresponds to the default constructed (invalid) vertex
  // accessor
  (void)level;
  (void)index;
  Assert((level == -2) && (index == -2),
         ExcMessage(
           "This constructor can not be called for face iterators in 1d, "
           "except to default-construct iterator objects."));
}
 
 
 
template <int spacedim>
template <int structdim2, int dim2, int spacedim2>
inline TriaAccessor<0, 1, spacedim>::TriaAccessor(
  const TriaAccessor<structdim2, dim2, spacedim2> &)
  : tria(nullptr)
  , vertex_kind(interior_vertex)
  , global_vertex_index(numbers::invalid_unsigned_int)
{
  Assert(false, ExcImpossibleInDim(0));
}
 
 
 
template <int spacedim>
template <int structdim2, int dim2, int spacedim2>
inline TriaAccessor<0, 1, spacedim>::TriaAccessor(
  const InvalidAccessor<structdim2, dim2, spacedim2> &)
  : tria(nullptr)
  , vertex_kind(interior_vertex)
  , global_vertex_index(numbers::invalid_unsigned_int)
{
  Assert(false, ExcImpossibleInDim(0));
}
 
 
 
template <int spacedim>
inline void
TriaAccessor<0, 1, spacedim>::copy_from(const TriaAccessor &t)
{
  tria                = t.tria;
  vertex_kind         = t.vertex_kind;
  global_vertex_index = t.global_vertex_index;
}
 
 
 
template <int spacedim>
inline void
TriaAccessor<0, 1, spacedim>::copy_from(
  const TriaAccessorBase<0, 1, spacedim> &)
{
  // We cannot convert from TriaAccessorBase to
  // TriaAccessor<0,1,spacedim> because the latter is not derived from
  // the former. We should never get here.
  DEAL_II_ASSERT_UNREACHABLE();
}
 
 
 
template <int spacedim>
inline bool
TriaAccessor<0, 1, spacedim>::operator<(
  const TriaAccessor<0, 1, spacedim> &other) const
{
  Assert(tria == other.tria, TriaAccessorExceptions::ExcCantCompareIterators());
 
  return (global_vertex_index < other.global_vertex_index);
}
 
 
 
template <int spacedim>
inline IteratorState::IteratorStates
TriaAccessor<0, 1, spacedim>::state()
{
  return IteratorState::valid;
}
 
 
template <int spacedim>
inline int
TriaAccessor<0, 1, spacedim>::level()
{
  return 0;
}
 
 
 
template <int spacedim>
inline int
TriaAccessor<0, 1, spacedim>::index() const
{
  return global_vertex_index;
}
 
 
 
template <int spacedim>
inline const Triangulation<1, spacedim> &
TriaAccessor<0, 1, spacedim>::get_triangulation() const
{
  return *tria;
}
 
 
 
template <int spacedim>
inline void
TriaAccessor<0, 1, spacedim>::operator++() const
{
  DEAL_II_NOT_IMPLEMENTED();
}
 
 
template <int spacedim>
inline void
TriaAccessor<0, 1, spacedim>::operator--() const
{
  DEAL_II_NOT_IMPLEMENTED();
}
 
 
 
template <int spacedim>
inline bool
TriaAccessor<0, 1, spacedim>::operator==(const TriaAccessor &t) const
{
  const bool result =
    ((tria == t.tria) && (global_vertex_index == t.global_vertex_index));
  // if we point to the same vertex, make sure we know the same about it
  if (result == true)
    Assert(vertex_kind == t.vertex_kind, ExcInternalError());
 
  return result;
}
 
 
 
template <int spacedim>
inline bool
TriaAccessor<0, 1, spacedim>::operator!=(const TriaAccessor &t) const
{
  return !(*this == t);
}
 
 
 
template <int spacedim>
inline unsigned int
TriaAccessor<0, 1, spacedim>::vertex_index(const unsigned int i) const
{
  AssertIndexRange(i, 1);
  (void)i;
  return global_vertex_index;
}
 
 
 
template <int spacedim>
inline Point<spacedim> &
TriaAccessor<0, 1, spacedim>::vertex(const unsigned int i) const
{
  AssertIndexRange(i, 1);
  (void)i;
  return const_cast<Point<spacedim> &>(
    this->tria->vertices[global_vertex_index]);
}
 
 
 
template <int spacedim>
inline Point<spacedim>
TriaAccessor<0, 1, spacedim>::center() const
{
  return this->tria->vertices[global_vertex_index];
}
 
 
 
template <int spacedim>
inline typename ::internal::TriangulationImplementation::
  Iterators<1, spacedim>::line_iterator
  TriaAccessor<0, 1, spacedim>::line(const unsigned int)
{
  return {};
}
 
 
template <int spacedim>
inline unsigned int
TriaAccessor<0, 1, spacedim>::line_index(const unsigned int)
{
  Assert(false, ExcImpossibleInDim(0));
  return numbers::invalid_unsigned_int;
}
 
 
template <int spacedim>
inline typename ::internal::TriangulationImplementation::
  Iterators<1, spacedim>::quad_iterator
  TriaAccessor<0, 1, spacedim>::quad(const unsigned int)
{
  return {};
}
 
 
 
template <int spacedim>
inline unsigned int
TriaAccessor<0, 1, spacedim>::quad_index(const unsigned int)
{
  Assert(false, ExcImpossibleInDim(0));
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int spacedim>
inline double
TriaAccessor<0, 1, spacedim>::measure()
{
  return 1.0;
}
 
 
 
template <int spacedim>
inline bool
TriaAccessor<0, 1, spacedim>::at_boundary() const
{
  return vertex_kind != interior_vertex;
}
 
 
template <int spacedim>
inline types::boundary_id
TriaAccessor<0, 1, spacedim>::boundary_id() const
{
  switch (vertex_kind)
    {
      case left_vertex:
      case right_vertex:
        {
          Assert(tria->vertex_to_boundary_id_map_1d->find(
                   this->vertex_index()) !=
                   tria->vertex_to_boundary_id_map_1d->end(),
                 ExcInternalError());
 
          return (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()];
        }
 
      default:
        return numbers::internal_face_boundary_id;
    }
}
 
 
 
template <int spacedim>
inline const Manifold<1, spacedim> &
TriaAccessor<0, 1, spacedim>::get_manifold() const
{
  return this->tria->get_manifold(this->manifold_id());
}
 
 
 
template <int spacedim>
inline types::manifold_id
TriaAccessor<0, 1, spacedim>::manifold_id() const
{
  if (tria->vertex_to_manifold_id_map_1d->find(this->vertex_index()) !=
      tria->vertex_to_manifold_id_map_1d->end())
    return (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()];
  else
    return numbers::flat_manifold_id;
}
 
 
template <int spacedim>
inline types::geometric_orientation
TriaAccessor<0, 1, spacedim>::combined_face_orientation(
  const unsigned int /*face*/)
{
  return numbers::reverse_line_orientation;
}
 
 
template <int spacedim>
inline bool
TriaAccessor<0, 1, spacedim>::face_orientation(const unsigned int /*face*/)
{
  return false;
}
 
 
 
template <int spacedim>
inline bool
TriaAccessor<0, 1, spacedim>::face_flip(const unsigned int /*face*/)
{
  return false;
}
 
 
 
template <int spacedim>
inline bool
TriaAccessor<0, 1, spacedim>::face_rotation(const unsigned int /*face*/)
{
  return false;
}
 
 
 
template <int spacedim>
inline types::geometric_orientation
TriaAccessor<0, 1, spacedim>::line_orientation(const unsigned int /*line*/)
{
  return numbers::reverse_line_orientation;
}
 
 
 
template <int spacedim>
inline bool
TriaAccessor<0, 1, spacedim>::has_children()
{
  return false;
}
 
 
 
template <int spacedim>
inline unsigned int
TriaAccessor<0, 1, spacedim>::n_children()
{
  return 0;
}
 
 
 
template <int spacedim>
inline unsigned int
TriaAccessor<0, 1, spacedim>::n_active_descendants()
{
  return 0;
}
 
 
 
template <int spacedim>
inline unsigned int
TriaAccessor<0, 1, spacedim>::max_refinement_depth()
{
  return 0;
}
 
 
 
template <int spacedim>
inline unsigned int
TriaAccessor<0, 1, spacedim>::child_iterator_to_index(
  const TriaIterator<TriaAccessor<0, 1, spacedim>> &)
{
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int spacedim>
inline TriaIterator<TriaAccessor<0, 1, spacedim>>
TriaAccessor<0, 1, spacedim>::child(const unsigned int)
{
  return TriaIterator<TriaAccessor<0, 1, spacedim>>();
}
 
 
template <int spacedim>
inline TriaIterator<TriaAccessor<0, 1, spacedim>>
TriaAccessor<0, 1, spacedim>::isotropic_child(const unsigned int)
{
  return TriaIterator<TriaAccessor<0, 1, spacedim>>();
}
 
 
template <int spacedim>
inline RefinementCase<0>
TriaAccessor<0, 1, spacedim>::refinement_case()
{
  return {RefinementPossibilities<0>::no_refinement};
}
 
template <int spacedim>
inline int
TriaAccessor<0, 1, spacedim>::child_index(const unsigned int)
{
  return -1;
}
 
 
template <int spacedim>
inline int
TriaAccessor<0, 1, spacedim>::isotropic_child_index(const unsigned int)
{
  return -1;
}
 
 
 
template <int spacedim>
inline void
TriaAccessor<0, 1, spacedim>::set_boundary_id(const types::boundary_id b) const
{
  Assert(tria->vertex_to_boundary_id_map_1d->find(this->vertex_index()) !=
           tria->vertex_to_boundary_id_map_1d->end(),
         ExcMessage("You can't set the boundary_id of a face of a cell that is "
                    "not actually at the boundary."));
 
  (*tria->vertex_to_boundary_id_map_1d)[this->vertex_index()] = b;
}
 
 
 
template <int spacedim>
inline void
TriaAccessor<0, 1, spacedim>::set_manifold_id(const types::manifold_id b)
{
  (*tria->vertex_to_manifold_id_map_1d)[this->vertex_index()] = b;
}
 
 
 
template <int spacedim>
inline void
TriaAccessor<0, 1, spacedim>::set_all_boundary_ids(
  const types::boundary_id b) const
{
  set_boundary_id(b);
}
 
 
 
template <int spacedim>
inline void
TriaAccessor<0, 1, spacedim>::set_all_manifold_ids(const types::manifold_id b)
{
  set_manifold_id(b);
}
 
 
 
template <int spacedim>
inline bool
TriaAccessor<0, 1, spacedim>::used() const
{
  return tria->vertex_used(global_vertex_index);
}
 
 
 
template <int spacedim>
inline ReferenceCell<0>
TriaAccessor<0, 1, spacedim>::reference_cell() const
{
  return ReferenceCells::Vertex;
}
 
 
 
template <int spacedim>
unsigned int
TriaAccessor<0, 1, spacedim>::n_vertices() const
{
  return 1;
}
 
 
 
template <int spacedim>
unsigned int
TriaAccessor<0, 1, spacedim>::n_lines() const
{
  return 0;
}
 
 
 
template <int spacedim>
std_cxx20::ranges::iota_view<unsigned int, unsigned int>
TriaAccessor<0, 1, spacedim>::vertex_indices() const
{
  return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
                                                                  n_vertices());
}
 
 
 
template <int spacedim>
std_cxx20::ranges::iota_view<unsigned int, unsigned int>
TriaAccessor<0, 1, spacedim>::line_indices() const
{
  return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(0U,
                                                                  n_lines());
}
 
/*------------------ Functions: CellAccessor<dim,spacedim> ------------------*/
 
 
template <int dim, int spacedim>
inline CellAccessor<dim, spacedim>::CellAccessor(
  const Triangulation<dim, spacedim> *parent,
  const int                           level,
  const int                           index,
  const AccessorData                 *local_data)
  : TriaAccessor<dim, dim, spacedim>(parent, level, index, local_data)
{}
 
 
 
template <int dim, int spacedim>
inline CellAccessor<dim, spacedim>::CellAccessor(
  const TriaAccessor<dim, dim, spacedim> &cell_accessor)
  : TriaAccessor<dim, dim, spacedim>(
      static_cast<const TriaAccessor<dim, dim, spacedim> &>(cell_accessor))
{}
 
 
 
template <int dim, int spacedim>
inline TriaIterator<CellAccessor<dim, spacedim>>
CellAccessor<dim, spacedim>::child(const unsigned int i) const
{
  TriaIterator<CellAccessor<dim, spacedim>> q(this->tria,
                                              this->level() + 1,
                                              this->child_index(i));
 
  Assert((q.state() == IteratorState::past_the_end) || q->used(),
         ExcInternalError());
 
  return q;
}
 
 
 
template <int dim, int spacedim>
inline boost::container::small_vector<TriaIterator<CellAccessor<dim, spacedim>>,
                                      ReferenceCells::max_n_children<dim>()>
CellAccessor<dim, spacedim>::child_iterators() const
{
  boost::container::small_vector<TriaIterator<CellAccessor<dim, spacedim>>,
                                 ReferenceCells::max_n_children<dim>()>
    child_iterators(this->n_children());
 
  for (unsigned int i = 0; i < this->n_children(); ++i)
    child_iterators[i] = this->child(i);
 
  return child_iterators;
}
 
 
 
template <int dim, int spacedim>
inline TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>
CellAccessor<dim, spacedim>::face(const unsigned int i) const
{
  AssertIndexRange(i, this->n_faces());
  if constexpr (dim == 1)
    {
      using VertexKind = typename TriaAccessor<0, 1, spacedim>::VertexKind;
      VertexKind vertex_kind = VertexKind::interior_vertex;
      if (i == 0 && at_boundary(0))
        vertex_kind = VertexKind::left_vertex;
      if (i == 1 && at_boundary(1))
        vertex_kind = VertexKind::right_vertex;
      TriaAccessor<0, 1, spacedim> a(&this->get_triangulation(),
                                     vertex_kind,
                                     this->vertex_index(i));
      return TriaIterator<TriaAccessor<0, 1, spacedim>>(a);
    }
  else if constexpr (dim == 2)
    return this->line(i);
  else if constexpr (dim == 3)
    return this->quad(i);
  else
    {
      Assert(false, ExcNotImplemented());
      return {};
    }
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
CellAccessor<dim, spacedim>::face_iterator_to_index(
  const TriaIterator<TriaAccessor<dim - 1, dim, spacedim>> &face) const
{
  for (const unsigned int face_n : this->face_indices())
    if (this->face(face_n) == face)
      return face_n;
 
  Assert(false,
         ExcMessage("The given face is not a face of the current cell."));
  return numbers::invalid_unsigned_int;
}
 
 
 
template <int dim, int spacedim>
inline boost::container::small_vector<
  TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
#  ifndef _MSC_VER
  ReferenceCells::max_n_faces<dim>()
#  else
  GeometryInfo<dim>::faces_per_cell
#  endif
  >
CellAccessor<dim, spacedim>::face_iterators() const
{
  boost::container::small_vector<
    TriaIterator<TriaAccessor<dim - 1, dim, spacedim>>,
#  ifndef _MSC_VER
    ReferenceCells::max_n_faces<dim>()
#  else
    GeometryInfo<dim>::faces_per_cell
#  endif
    >
    face_iterators(this->n_faces());
 
  for (const unsigned int i : this->face_indices())
    face_iterators[i] = this->face(i);
 
  return face_iterators;
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
CellAccessor<dim, spacedim>::face_index(const unsigned int i) const
{
  switch (dim)
    {
      case 1:
        return this->vertex_index(i);
 
      case 2:
        return this->line_index(i);
 
      case 3:
        return this->quad_index(i);
 
      default:
        return numbers::invalid_unsigned_int;
    }
}
 
 
 
template <int dim, int spacedim>
inline int
CellAccessor<dim, spacedim>::neighbor_index(const unsigned int face_no) const
{
  AssertIndexRange(face_no, this->n_faces());
  return this->tria->levels[this->level()]->neighbor_index(this->present_index,
                                                           face_no);
}
 
 
 
template <int dim, int spacedim>
inline int
CellAccessor<dim, spacedim>::neighbor_level(const unsigned int face_no) const
{
  AssertIndexRange(face_no, this->n_faces());
  return this->tria->levels[this->level()]->neighbor_level(this->present_index,
                                                           face_no);
}
 
 
 
template <int dim, int spacedim>
inline RefinementCase<dim>
CellAccessor<dim, spacedim>::refine_flag_set() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  // cells flagged for refinement must be active (the @p set_refine_flag
  // function checks this, but activity may change when refinement is executed
  // and for some reason the refine flag is not cleared).
  Assert(
    this->is_active() ||
      !this->tria->levels[this->level()]->refine_flags[this->present_index],
    ExcRefineCellNotActive());
  return RefinementCase<dim>(
    this->tria->levels[this->level()]->refine_flags[this->present_index]);
}
 
 
 
template <int dim, int spacedim>
inline void
CellAccessor<dim, spacedim>::set_refine_flag(
  const RefinementCase<dim> refinement_case) const
{
  Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
  Assert(!coarsen_flag_set(), ExcCellFlaggedForCoarsening());
 
  this->tria->levels[this->level()]->refine_flags[this->present_index] =
    refinement_case;
}
 
 
 
template <int dim, int spacedim>
inline void
CellAccessor<dim, spacedim>::clear_refine_flag() const
{
  Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
  this->tria->levels[this->level()]->refine_flags[this->present_index] =
    RefinementCase<dim>::no_refinement;
}
 
 
template <int dim, int spacedim>
inline std::uint8_t
CellAccessor<dim, spacedim>::refine_choice() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  if (this->tria->levels[this->level()]->refine_choice.size() == 0)
    return 0U;
  return this->tria->levels[this->level()]->refine_choice[this->present_index];
}
 
 
template <int dim, int spacedim>
inline void
CellAccessor<dim, spacedim>::set_refine_choice(
  const std::uint8_t refinement_choice) const
{
  Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
  if (this->tria->levels[this->level()]->refine_choice.size() != 0)
    this->tria->levels[this->level()]->refine_choice[this->present_index] =
      refinement_choice;
}
 
 
template <int dim, int spacedim>
inline void
CellAccessor<dim, spacedim>::clear_refine_choice() const
{
  Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
  if (this->tria->levels[this->level()]->refine_choice.size() != 0)
    this->tria->levels[this->level()]->refine_choice[this->present_index] =
      static_cast<char>(IsotropicRefinementChoice::isotropic_refinement);
}
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::flag_for_face_refinement(
  const unsigned int             face_no,
  const RefinementCase<dim - 1> &face_refinement_case) const
{
  Assert(dim > 1, ExcImpossibleInDim(dim));
  AssertIndexRange(face_no, this->n_faces());
  AssertIndexRange(face_refinement_case,
                   RefinementCase<dim>::isotropic_refinement + 1);
 
  // the new refinement case is a combination of the minimum required one for
  // the given face refinement and the already existing flagged refinement case
  RefinementCase<dim> old_ref_case = refine_flag_set();
  RefinementCase<dim> new_ref_case =
    (old_ref_case |
     GeometryInfo<dim>::min_cell_refinement_case_for_face_refinement(
       face_refinement_case,
       face_no,
       this->face_orientation(face_no),
       this->face_flip(face_no),
       this->face_rotation(face_no)));
  set_refine_flag(new_ref_case);
  // return, whether we had to change the refinement flag
  return new_ref_case != old_ref_case;
}
 
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::flag_for_line_refinement(
  const unsigned int line_no) const
{
  Assert(dim > 1, ExcImpossibleInDim(dim));
  AssertIndexRange(line_no, this->n_lines());
 
  // the new refinement case is a combination of the minimum required one for
  // the given line refinement and the already existing flagged refinement case
  RefinementCase<dim>
    old_ref_case = refine_flag_set(),
    new_ref_case =
      old_ref_case |
      GeometryInfo<dim>::min_cell_refinement_case_for_line_refinement(line_no);
  set_refine_flag(new_ref_case);
  // return, whether we had to change the refinement flag
  return new_ref_case != old_ref_case;
}
 
 
 
template <int dim, int spacedim>
inline internal::SubfaceCase<dim>
CellAccessor<dim, spacedim>::subface_case(const unsigned int face_no) const
{
  Assert(is_active(), TriaAccessorExceptions::ExcCellNotActive());
  AssertIndexRange(face_no, this->n_faces());
 
  if constexpr (dim == 1)
    return internal::SubfaceCase<1>::case_none;
  else if constexpr (dim == 2)
    return ((face(face_no)->has_children()) ?
              internal::SubfaceCase<2>::case_x :
              internal::SubfaceCase<2>::case_none);
  else if constexpr (dim == 3)
    {
      switch (static_cast<std::uint8_t>(face(face_no)->refinement_case()))
        {
          case RefinementCase<3>::no_refinement:
            return internal::SubfaceCase<3>::case_none;
          case RefinementCase<3>::cut_x:
            if (face(face_no)->child(0)->has_children())
              {
                Assert(face(face_no)->child(0)->refinement_case() ==
                         RefinementCase<2>::cut_y,
                       ExcInternalError());
                if (face(face_no)->child(1)->has_children())
                  {
                    Assert(face(face_no)->child(1)->refinement_case() ==
                             RefinementCase<2>::cut_y,
                           ExcInternalError());
                    return internal::SubfaceCase<3>::case_x1y2y;
                  }
                else
                  return internal::SubfaceCase<3>::case_x1y;
              }
            else
              {
                if (face(face_no)->child(1)->has_children())
                  {
                    Assert(face(face_no)->child(1)->refinement_case() ==
                             RefinementCase<2>::cut_y,
                           ExcInternalError());
                    return internal::SubfaceCase<3>::case_x2y;
                  }
                else
                  return internal::SubfaceCase<3>::case_x;
              }
          case RefinementCase<3>::cut_y:
            if (face(face_no)->child(0)->has_children())
              {
                Assert(face(face_no)->child(0)->refinement_case() ==
                         RefinementCase<2>::cut_x,
                       ExcInternalError());
                if (face(face_no)->child(1)->has_children())
                  {
                    Assert(face(face_no)->child(1)->refinement_case() ==
                             RefinementCase<2>::cut_x,
                           ExcInternalError());
                    return internal::SubfaceCase<3>::case_y1x2x;
                  }
                else
                  return internal::SubfaceCase<3>::case_y1x;
              }
            else
              {
                if (face(face_no)->child(1)->has_children())
                  {
                    Assert(face(face_no)->child(1)->refinement_case() ==
                             RefinementCase<2>::cut_x,
                           ExcInternalError());
                    return internal::SubfaceCase<3>::case_y2x;
                  }
                else
                  return internal::SubfaceCase<3>::case_y;
              }
          case RefinementCase<3>::cut_xy:
            return internal::SubfaceCase<3>::case_xy;
          default:
            DEAL_II_ASSERT_UNREACHABLE();
        }
    }
 
  // we should never get here
  DEAL_II_ASSERT_UNREACHABLE();
  return internal::SubfaceCase<dim>::case_none;
}
 
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::coarsen_flag_set() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  // cells flagged for coarsening must be active (the @p set_refine_flag
  // function checks this, but activity may change when refinement is executed
  // and for some reason the refine flag is not cleared).
  Assert(
    this->is_active() ||
      !this->tria->levels[this->level()]->coarsen_flags[this->present_index],
    ExcRefineCellNotActive());
  return this->tria->levels[this->level()]->coarsen_flags[this->present_index];
}
 
 
 
template <int dim, int spacedim>
inline void
CellAccessor<dim, spacedim>::set_coarsen_flag() const
{
  Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
  Assert(!refine_flag_set(), ExcCellFlaggedForRefinement());
 
  this->tria->levels[this->level()]->coarsen_flags[this->present_index] = true;
}
 
 
 
template <int dim, int spacedim>
inline void
CellAccessor<dim, spacedim>::clear_coarsen_flag() const
{
  Assert(this->used() && this->is_active(), ExcRefineCellNotActive());
  this->tria->levels[this->level()]->coarsen_flags[this->present_index] = false;
}
 
 
 
template <int dim, int spacedim>
inline TriaIterator<CellAccessor<dim, spacedim>>
CellAccessor<dim, spacedim>::neighbor(const unsigned int face_no) const
{
  TriaIterator<CellAccessor<dim, spacedim>> q(this->tria,
                                              neighbor_level(face_no),
                                              neighbor_index(face_no));
 
  Assert((q.state() == IteratorState::past_the_end) || q->used(),
         ExcInternalError());
 
  return q;
}
 
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::is_active() const
{
  return !this->has_children();
}
 
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::is_locally_owned() const
{
  Assert(this->is_active(),
         ExcMessage("is_locally_owned() can only be called on active cells!"));
#  ifndef DEAL_II_WITH_MPI
  return true;
#  else
 
  // Serial triangulations report invalid_subdomain_id as their locally owned
  // subdomain, so the first condition checks whether we have a serial
  // triangulation, in which case all cells are locally owned. The second
  // condition compares the subdomain id in the parallel case.
  const types::subdomain_id locally_owned_subdomain =
    this->tria->locally_owned_subdomain();
  return (locally_owned_subdomain == numbers::invalid_subdomain_id ||
          this->subdomain_id() == locally_owned_subdomain);
 
#  endif
}
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::is_locally_owned_on_level() const
{
#  ifndef DEAL_II_WITH_MPI
  return true;
#  else
 
  // Serial triangulations report invalid_subdomain_id as their locally owned
  // subdomain, so the first condition checks whether we have a serial
  // triangulation, in which case all cells are locally owned. The second
  // condition compares the subdomain id in the parallel case.
  const types::subdomain_id locally_owned_subdomain =
    this->tria->locally_owned_subdomain();
  return (locally_owned_subdomain == numbers::invalid_subdomain_id ||
          this->level_subdomain_id() == locally_owned_subdomain);
 
#  endif
}
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::is_ghost() const
{
  Assert(this->is_active(),
         ExcMessage("is_ghost() can only be called on active cells!"));
  if (this->has_children())
    return false;
 
#  ifndef DEAL_II_WITH_MPI
  return false;
#  else
 
  // Serial triangulations report invalid_subdomain_id as their locally owned
  // subdomain, so the first condition rules out that case as all cells to a
  // serial triangulation are locally owned and none is ghosted. The second
  // and third conditions check whether the cell's subdomain is not the
  // locally owned one and not artificial.
  const types::subdomain_id locally_owned_subdomain =
    this->tria->locally_owned_subdomain();
  const types::subdomain_id subdomain_id = this->subdomain_id();
  return (locally_owned_subdomain != numbers::invalid_subdomain_id &&
          subdomain_id != locally_owned_subdomain &&
          subdomain_id != numbers::artificial_subdomain_id);
 
#  endif
}
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::is_ghost_on_level() const
{
#  ifndef DEAL_II_WITH_MPI
  return false;
#  else
 
  // Serial triangulations report invalid_subdomain_id as their locally owned
  // subdomain, so the first condition checks whether we have a serial
  // triangulation, in which case all cells are locally owned. The second
  // condition compares the subdomain id in the parallel case.
  const types::subdomain_id locally_owned_subdomain =
    this->tria->locally_owned_subdomain();
  const types::subdomain_id subdomain_id = this->level_subdomain_id();
  return (locally_owned_subdomain != numbers::invalid_subdomain_id &&
          subdomain_id != locally_owned_subdomain &&
          subdomain_id != numbers::artificial_subdomain_id);
 
#  endif
}
 
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::is_artificial() const
{
  Assert(this->is_active(),
         ExcMessage("is_artificial() can only be called on active cells!"));
#  ifndef DEAL_II_WITH_MPI
  return false;
#  else
 
  // Serial triangulations report invalid_subdomain_id as their locally owned
  // subdomain, so the first condition rules out that case as all cells to a
  // serial triangulation are locally owned and none is artificial.
  return (this->tria->locally_owned_subdomain() !=
            numbers::invalid_subdomain_id &&
          this->subdomain_id() == numbers::artificial_subdomain_id);
 
#  endif
}
 
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::is_artificial_on_level() const
{
#  ifndef DEAL_II_WITH_MPI
  return false;
#  else
  return (this->tria->locally_owned_subdomain() !=
            numbers::invalid_subdomain_id &&
          this->level_subdomain_id() == numbers::artificial_subdomain_id);
#  endif
}
 
 
 
template <int dim, int spacedim>
inline void
CellAccessor<dim, spacedim>::set_neighbor(
  const unsigned int                               face_no,
  const TriaIterator<CellAccessor<dim, spacedim>> &neighbor) const
{
  AssertIndexRange(face_no, this->n_faces());
 
  Assert(neighbor.state() != IteratorState::invalid,
         ExcMessage("The provided neighbor must be either a valid iterator or "
                    "an end iterator (to indicate a boundary face)"));
  auto &level = this->tria->levels[this->level()];
  level->set_neighbor(this->present_index,
                      face_no,
                      neighbor->level(),
                      neighbor->index());
}
 
 
 
template <int dim, int spacedim>
inline types::material_id
CellAccessor<dim, spacedim>::material_id() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  return this->tria->levels[this->level()]
    ->cells.boundary_or_material_id[this->present_index]
    .material_id;
}
 
 
 
template <int dim, int spacedim>
inline types::subdomain_id
CellAccessor<dim, spacedim>::subdomain_id() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  Assert(this->is_active(),
         ExcMessage("subdomain_id() can only be called on active cells!"));
  return this->tria->levels[this->level()]->subdomain_ids[this->present_index];
}
 
 
 
template <int dim, int spacedim>
inline types::subdomain_id
CellAccessor<dim, spacedim>::level_subdomain_id() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  return this->tria->levels[this->level()]
    ->level_subdomain_ids[this->present_index];
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
CellAccessor<dim, spacedim>::neighbor_face_no(const unsigned int neighbor) const
{
  const unsigned int n2 = neighbor_of_neighbor_internal(neighbor);
  if (n2 != numbers::invalid_unsigned_int)
    // return this value as the neighbor is not coarser
    return n2;
  else
    // the neighbor is coarser
    return neighbor_of_coarser_neighbor(neighbor).first;
}
 
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::is_level_cell()
{
  return false;
}
 
 
 
template <int dim, int spacedim>
inline unsigned int
CellAccessor<dim, spacedim>::active_cell_index() const
{
  Assert(this->is_active(), TriaAccessorExceptions::ExcCellNotActive());
  return this->tria->levels[this->level()]
    ->active_cell_indices[this->present_index];
}
 
 
 
template <int dim, int spacedim>
inline types::global_cell_index
CellAccessor<dim, spacedim>::global_active_cell_index() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  Assert(this->is_active(),
         ExcMessage(
           "global_active_cell_index() can only be called on active cells!"));
 
  return this->tria->levels[this->level()]
    ->global_active_cell_indices[this->present_index];
}
 
 
 
template <int dim, int spacedim>
inline types::global_cell_index
CellAccessor<dim, spacedim>::global_level_cell_index() const
{
  return this->tria->levels[this->level()]
    ->global_level_cell_indices[this->present_index];
}
 
 
 
template <int dim, int spacedim>
inline bool
CellAccessor<dim, spacedim>::direction_flag() const
{
  Assert(this->used(), TriaAccessorExceptions::ExcCellNotUsed());
  if constexpr (dim == spacedim)
    return true;
  else if constexpr (dim == spacedim - 1)
    return this->tria->levels[this->level()]
      ->direction_flags[this->present_index];
  else
    {
      Assert(false,
             ExcMessage("This function cannot be called if dim<spacedim-1."));
      return true;
    }
}
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
 
#endif