]> https://gitweb.dealii.org/ - dealii.git/commitdiff
Don't use a .templates.h file if we unconditionally include it anyway.
authorWolfgang Bangerth <bangerth@colostate.edu>
Wed, 15 Jan 2025 03:39:14 +0000 (20:39 -0700)
committerWolfgang Bangerth <bangerth@colostate.edu>
Wed, 15 Jan 2025 14:53:35 +0000 (07:53 -0700)
include/deal.II/dofs/dof_accessor.h
include/deal.II/dofs/dof_accessor.templates.h

index 6921a50e1cc706239ef174ee89749455b1a19711..1584668420168f78df93267c1ad5bf8ebebd5b80 100644 (file)
 
 #include <deal.II/base/config.h>
 
+#include <deal.II/base/types.h>
+
+#include <deal.II/dofs/dof_faces.h>
 #include <deal.II/dofs/dof_handler.h>
 #include <deal.II/dofs/dof_iterator_selector.h>
+#include <deal.II/dofs/dof_levels.h>
 
 #include <deal.II/grid/tria_accessor.h>
 
+#include <deal.II/lac/affine_constraints.h>
 #include <deal.II/lac/read_vector.h>
+#include <deal.II/lac/read_write_vector.h>
 
 #include <boost/container/small_vector.hpp>
 
+#include <limits>
 #include <set>
+#include <type_traits>
 #include <vector>
 
+
 DEAL_II_NAMESPACE_OPEN
 
 // Forward declarations
@@ -73,11 +82,6 @@ namespace internal
 } // namespace internal
 #endif
 
-// note: the file dof_accessor.templates.h is included at the end of
-// this file.  this includes a lot of templates and thus makes
-// compilation slower, but at the same time allows for more aggressive
-// inlining and thus faster code.
-
 
 namespace internal
 {
@@ -2139,10 +2143,2947 @@ DoFInvalidAccessor<structdim, dim, spacedim>::DoFInvalidAccessor(
 
 
 
-DEAL_II_NAMESPACE_CLOSE
+/*------------------------- Functions: DoFAccessor ---------------------------*/
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline DoFAccessor<structdim, dim, spacedim, level_dof_access>::DoFAccessor()
+{
+  Assert(false, ExcInvalidObject());
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline DoFAccessor<structdim, dim, spacedim, level_dof_access>::DoFAccessor(
+  const Triangulation<dim, spacedim> *tria,
+  const int                           level,
+  const int                           index,
+  const DoFHandler<dim, spacedim>    *dof_handler)
+  : dealii::internal::DoFAccessorImplementation::
+      Inheritance<structdim, dim, spacedim>::BaseClass(tria, level, index)
+  , dof_handler(const_cast<DoFHandler<dim, spacedim> *>(dof_handler))
+{
+  Assert(
+    tria == nullptr || &dof_handler->get_triangulation() == tria,
+    ExcMessage(
+      "You can't create a DoF accessor in which the DoFHandler object "
+      "uses a different triangulation than the one you pass as argument."));
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+template <int structdim2, int dim2, int spacedim2>
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::DoFAccessor(
+  const InvalidAccessor<structdim2, dim2, spacedim2> &)
+{
+  Assert(false, ExcInvalidObject());
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
+inline DoFAccessor<structdim, dim, spacedim, level_dof_access>::DoFAccessor(
+  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &other)
+  : BaseClass(other)
+  , dof_handler(nullptr)
+{
+  Assert(false,
+         ExcMessage(
+           "You are trying to assign iterators that are incompatible. "
+           "The reason for incompatibility is that they refer to objects of "
+           "different dimensionality (e.g., assigning a line iterator "
+           "to a quad iterator)."));
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+template <bool level_dof_access2>
+inline DoFAccessor<structdim, dim, spacedim, level_dof_access>::DoFAccessor(
+  const DoFAccessor<structdim, dim, spacedim, level_dof_access2> &other)
+  : BaseClass(other)
+  , dof_handler(const_cast<DoFHandler<dim, spacedim> *>(other.dof_handler))
+{}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline void
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::set_dof_handler(
+  DoFHandler<dim, spacedim> *dh)
+{
+  Assert(dh != nullptr, ExcInvalidObject());
+  this->dof_handler = dh;
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline const DoFHandler<dim, spacedim> &
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::get_dof_handler() const
+{
+  Assert(this->dof_handler != nullptr, ExcInvalidObject());
+  return *this->dof_handler;
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline void
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::copy_from(
+  const TriaAccessorBase<structdim, dim, spacedim> &da)
+{
+  Assert(this->dof_handler != nullptr, ExcInvalidObject());
+  BaseClass::copy_from(da);
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+template <bool level_dof_access2>
+inline void
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::copy_from(
+  const DoFAccessor<structdim, dim, spacedim, level_dof_access2> &a)
+{
+  BaseClass::copy_from(a);
+  this->dof_handler = a.dof_handler;
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
+inline bool
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::operator==(
+  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &a) const
+{
+  Assert(structdim == structdim2, ExcCantCompareIterators());
+  Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
+  return (BaseClass::operator==(a));
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
+inline bool
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::operator!=(
+  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &a) const
+{
+  Assert(structdim == structdim2, ExcCantCompareIterators());
+  Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
+  return (BaseClass::operator!=(a));
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline TriaIterator<DoFAccessor<structdim, dim, spacedim, level_dof_access>>
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::child(
+  const unsigned int i) const
+{
+  Assert(static_cast<unsigned int>(this->present_level) <
+           this->dof_handler->object_dof_indices.size(),
+         ExcMessage("DoFHandler not initialized"));
+
+  TriaIterator<TriaAccessor<structdim, dim, spacedim>> t =
+    TriaAccessor<structdim, dim, spacedim>::child(i);
+
+  TriaIterator<DoFAccessor<structdim, dim, spacedim, level_dof_access>> q(
+    *t, this->dof_handler);
+  return q;
+}
+
+
+namespace internal
+{
+  namespace DoFAccessorImplementation
+  {
+    /**
+     * Convert an FE index that might contain the right value but also
+     * numbers::invalid_fe_index to a right value if needed/possible.
+     */
+    template <int structdim, int dim, int spacedim, bool level_dof_access>
+    types::fe_index
+    get_fe_index_or_default(
+      const DoFAccessor<structdim, dim, spacedim, level_dof_access> &cell,
+      const types::fe_index                                          fe_index)
+    {
+      if (cell.get_dof_handler().has_hp_capabilities() == false)
+        {
+          // No hp enabled, and the argument is at its default value -> we
+          // can translate to the default active fe index
+          Assert(
+            (fe_index == numbers::invalid_fe_index) ||
+              (fe_index == DoFHandler<dim, spacedim>::default_fe_index),
+            ExcMessage(
+              "It is not possible to specify a FE index if no hp support is used!"));
+
+          return DoFHandler<dim, spacedim>::default_fe_index;
+        }
+      else
+        {
+          // Otherwise: If anything other than the default is provided by
+          // the caller, then we should take just that. As an exception, if
+          // we are on a cell (rather than a face/edge/vertex), then we know
+          // that there is only one active fe index on this cell and we can
+          // use that:
+          if ((dim == structdim) && (fe_index == numbers::invalid_fe_index))
+            {
+              AssertDimension(cell.n_active_fe_indices(), 1);
+
+              return cell.nth_active_fe_index(0);
+            }
+
+          Assert((fe_index != numbers::invalid_fe_index),
+                 ExcMessage(
+                   "You need to specify a FE index if hp support is used!"));
+
+          return fe_index;
+        }
+    }
+
+    /**
+     * A class like the one with same name in tria.cc. See there for more
+     * information.
+     */
+    struct Implementation
+    {
+      /**
+       * In several applications of DoFAccessor::get_dof_values(), we want to
+       * extract some indices without having to allocate memory. We do this by
+       * setting a boost small_vector with 27 elements on the stack, and only
+       * allocate when we exceed 27. The number 27 is heuristic and allows up
+       * to quadratic shape functions on scalar problems in 3d, or linear
+       * shape functions on systems (elasticity).
+       */
+      using dof_index_vector_type =
+        boost::container::small_vector<dealii::types::global_dof_index, 27>;
+
+      /**
+       * Process the @p local_index-th degree of freedom corresponding to the
+       * finite element specified by @p fe_index on the vertex with global
+       * number @p vertex_index to @p global_index.
+       *
+       * The template argument `structdim` indicates the
+       * dimensionality of the object on which we seek to know the DoF
+       * index. For example, if `structdim==0`, then we are looking to
+       * get a DoF index on a vertex of the indicated cell.
+       */
+      template <int dim,
+                int spacedim,
+                int structdim,
+                typename GlobalIndexType,
+                typename DoFPProcessor>
+      static void
+      process_dof_index(const DoFHandler<dim, spacedim> &dof_handler,
+                        const unsigned int               obj_level,
+                        const unsigned int               obj_index,
+                        const types::fe_index            fe_index,
+                        const unsigned int               local_index,
+                        const std::integral_constant<int, structdim> &,
+                        GlobalIndexType     &global_index,
+                        const DoFPProcessor &process)
+      {
+        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
+
+        // 1) no hp used -> fe_index == 0
+        if (dof_handler.hp_capability_enabled == false)
+          {
+            AssertDimension(fe_index,
+                            (DoFHandler<dim, spacedim>::default_fe_index));
+
+            process(
+              dof_handler.object_dof_indices
+                [obj_level][structdim]
+                [dof_handler.object_dof_ptr[obj_level][structdim][obj_index] +
+                 local_index],
+              global_index);
+
+            return;
+          }
+
+        // 2) cell and hp is used -> there is only one fe_index
+        if (structdim == dim)
+          {
+            process(
+              dof_handler.object_dof_indices
+                [obj_level][structdim]
+                [dof_handler.object_dof_ptr[obj_level][structdim][obj_index] +
+                 local_index],
+              global_index);
+            return;
+          }
+
+        // 3) general entity and hp is used
+        AssertIndexRange(obj_level, dof_handler.object_dof_indices.size());
+        AssertIndexRange(structdim,
+                         dof_handler.object_dof_indices[obj_level].size());
+
+        Assert(dof_handler.hp_capability_enabled, ExcInternalError());
+
+        AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
+        AssertIndexRange(obj_index,
+                         dof_handler.hp_object_fe_ptr[structdim].size());
+
+        const auto ptr =
+          std::find(dof_handler.hp_object_fe_indices[structdim].begin() +
+                      dof_handler.hp_object_fe_ptr[structdim][obj_index],
+                    dof_handler.hp_object_fe_indices[structdim].begin() +
+                      dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
+                    fe_index);
+
+        Assert(ptr != dof_handler.hp_object_fe_indices[structdim].begin() +
+                        dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
+               ExcMessage(
+                 "You are requesting an active FE index that is not assigned "
+                 "to any of the cells connected to this entity."));
+
+        const types::fe_index fe_index_ =
+          std::distance(dof_handler.hp_object_fe_indices[structdim].begin() +
+                          dof_handler.hp_object_fe_ptr[structdim][obj_index],
+                        ptr);
+
+        AssertIndexRange(
+          dof_handler.hp_capability_enabled ?
+            (dof_handler.hp_object_fe_ptr[structdim][obj_index] + fe_index_) :
+            obj_index,
+          dof_handler.object_dof_ptr[obj_level][structdim].size());
+
+        AssertIndexRange(
+          dof_handler.object_dof_ptr
+              [obj_level][structdim]
+              [dof_handler.hp_capability_enabled ?
+                 (dof_handler.hp_object_fe_ptr[structdim][obj_index] +
+                  fe_index_) :
+                 obj_index] +
+            local_index,
+          dof_handler.object_dof_indices[obj_level][structdim].size());
+
+        process(dof_handler.object_dof_indices
+                  [obj_level][structdim]
+                  [dof_handler.object_dof_ptr
+                     [obj_level][structdim]
+                     [dof_handler.hp_capability_enabled ?
+                        (dof_handler.hp_object_fe_ptr[structdim][obj_index] +
+                         fe_index_) :
+                        obj_index] +
+                   local_index],
+                global_index);
+      }
+
+      /**
+       * Determine start index and number of dofs of object in global data
+       * structure.
+       *
+       * The template argument `structdim` indicates the
+       * dimensionality of the object on which we seek to know the DoF
+       * index. For example, if `structdim==0`, then we are looking to
+       * get a DoF index on a vertex of the indicated cell.
+       */
+      template <int dim, int spacedim, int structdim>
+      static std::pair<unsigned int, unsigned int>
+      process_object_range(const DoFHandler<dim, spacedim> &dof_handler,
+                           const unsigned int               obj_level,
+                           const unsigned int               obj_index,
+                           const types::fe_index            fe_index,
+                           const std::integral_constant<int, structdim> &)
+      {
+        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
+
+        // determine range of dofs in global data structure
+        // 1) cell
+        if (structdim == dim)
+          {
+            const unsigned int ptr_0 =
+              dof_handler.object_dof_ptr[obj_level][structdim][obj_index];
+            const unsigned int length =
+              dof_handler.get_fe(fe_index).template n_dofs_per_object<dim>(0);
+
+            return {ptr_0, length};
+          }
+
+        // 2) hp is not used -> fe_index == 0
+        if (dof_handler.hp_capability_enabled == false)
+          {
+            AssertDimension(fe_index,
+                            (DoFHandler<dim, spacedim>::default_fe_index));
+
+            const unsigned int ptr_0 =
+              dof_handler.object_dof_ptr[obj_level][structdim][obj_index];
+            const unsigned int length =
+              dof_handler.object_dof_ptr[obj_level][structdim][obj_index + 1] -
+              ptr_0;
+
+            return {ptr_0, length};
+          }
+
+        // 3) hp is used
+        AssertIndexRange(obj_level, dof_handler.object_dof_indices.size());
+        AssertIndexRange(structdim,
+                         dof_handler.object_dof_indices[obj_level].size());
+
+        AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
+        AssertIndexRange(obj_index,
+                         dof_handler.hp_object_fe_ptr[structdim].size());
+
+        const auto fe_index_local_ptr =
+          std::find(dof_handler.hp_object_fe_indices[structdim].begin() +
+                      dof_handler.hp_object_fe_ptr[structdim][obj_index],
+                    dof_handler.hp_object_fe_indices[structdim].begin() +
+                      dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
+                    fe_index);
+
+        Assert(fe_index_local_ptr !=
+                 dof_handler.hp_object_fe_indices[structdim].begin() +
+                   dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
+               ExcMessage(
+                 "You tried to call a function accessing DoF indices, but "
+                 "they appear not be available (yet) or inconsistent. "
+                 "Did you call distribute_dofs() first? Alternatively, if "
+                 "you are using different elements on different cells (i.e., "
+                 "you are using the hp capabilities of deal.II), did you "
+                 "change the active_fe_index of a cell since you last "
+                 "called distribute_dofs()?"));
+
+        const types::fe_index fe_index_local =
+          std::distance(dof_handler.hp_object_fe_indices[structdim].begin() +
+                          dof_handler.hp_object_fe_ptr[structdim][obj_index],
+                        fe_index_local_ptr);
+
+        AssertIndexRange(
+          dof_handler.hp_object_fe_ptr[structdim][obj_index] + fe_index_local,
+          dof_handler.object_dof_ptr[obj_level][structdim].size());
+
+        const unsigned int ptr_0 =
+          dof_handler
+            .object_dof_ptr[obj_level][structdim]
+                           [dof_handler.hp_object_fe_ptr[structdim][obj_index] +
+                            fe_index_local];
+        const unsigned int ptr_1 =
+          dof_handler
+            .object_dof_ptr[obj_level][structdim]
+                           [dof_handler.hp_object_fe_ptr[structdim][obj_index] +
+                            fe_index_local + 1];
+
+        return {ptr_0, ptr_1 - ptr_0};
+      }
+
+      template <int dim, int spacedim, int structdim, bool level_dof_access>
+      static std::pair<unsigned int, unsigned int>
+      process_object_range(
+        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
+                              accessor,
+        const types::fe_index fe_index)
+      {
+        return process_object_range(accessor.get_dof_handler(),
+                                    accessor.level(),
+                                    accessor.index(),
+                                    fe_index,
+                                    std::integral_constant<int, structdim>());
+      }
+
+      template <int dim, int spacedim, int structdim>
+      static std::pair<unsigned int, unsigned int>
+      process_object_range(dealii::DoFInvalidAccessor<structdim, dim, spacedim>,
+                           const unsigned int)
+      {
+        DEAL_II_ASSERT_UNREACHABLE();
+
+        return {0, 0};
+      }
+
+
+
+      /**
+       * Process all dofs of an object.
+       *
+       * The template argument `structdim` indicates the
+       * dimensionality of the object on which we seek to know the DoF
+       * index. For example, if `structdim==0`, then we are looking to
+       * get a DoF index on a vertex of the indicated cell.
+       */
+      template <int dim,
+                int spacedim,
+                int structdim,
+                typename DoFProcessor,
+                typename DoFMapping>
+      static DEAL_II_ALWAYS_INLINE void
+      process_object(const DoFHandler<dim, spacedim>              &dof_handler,
+                     const unsigned int                            obj_level,
+                     const unsigned int                            obj_index,
+                     const types::fe_index                         fe_index,
+                     const DoFMapping                             &mapping,
+                     const std::integral_constant<int, structdim> &dd,
+                     types::global_dof_index *&dof_indices_ptr,
+                     const DoFProcessor       &process)
+      {
+        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
+
+        // determine range of dofs in global data structure
+        const auto range =
+          process_object_range(dof_handler, obj_level, obj_index, fe_index, dd);
+        if (range.second == 0)
+          return;
+
+        std::vector<types::global_dof_index> &object_dof_indices =
+          dof_handler
+            .object_dof_indices[structdim < dim ? 0 : obj_level][structdim];
+        AssertIndexRange(range.first, object_dof_indices.size());
+        types::global_dof_index *DEAL_II_RESTRICT stored_indices =
+          object_dof_indices.data() + range.first;
+
+        // process dofs
+        for (unsigned int i = 0; i < range.second; ++i)
+          {
+            process(
+              stored_indices[(structdim == 0 || structdim == dim) ? i :
+                                                                    mapping(i)],
+              dof_indices_ptr);
+            if (dof_indices_ptr != nullptr)
+              ++dof_indices_ptr;
+          }
+      }
+
+
+
+      /**
+       * Set the @p local_index-th degree of freedom corresponding to the
+       * finite element specified by @p fe_index on the vertex with global
+       * number @p vertex_index to @p global_index.
+       *
+       * The template argument `structdim` indicates the
+       * dimensionality of the object on which we seek to know the DoF
+       * index. For example, if `structdim==0`, then we are looking to
+       * get a DoF index on a vertex of the indicated cell.
+       */
+      template <int dim, int spacedim, int structdim>
+      static void
+      set_dof_index(const DoFHandler<dim, spacedim>              &dof_handler,
+                    const unsigned int                            obj_level,
+                    const unsigned int                            obj_index,
+                    const types::fe_index                         fe_index,
+                    const unsigned int                            local_index,
+                    const std::integral_constant<int, structdim> &dd,
+                    const types::global_dof_index                 global_index)
+      {
+        process_dof_index(dof_handler,
+                          obj_level,
+                          obj_index,
+                          fe_index,
+                          local_index,
+                          dd,
+                          global_index,
+                          [](auto &ptr, const auto &value) { ptr = value; });
+      }
+
+
+      /**
+       * Get the @p local_index-th degree of freedom corresponding to the
+       * finite element specified by @p fe_index on the vertex with global
+       * number @p vertex_index to @p global_index.
+       *
+       * The template argument `structdim` indicates the
+       * dimensionality of the object on which we seek to know the DoF
+       * index. For example, if `structdim==0`, then we are looking to
+       * get a DoF index on a vertex of the indicated cell.
+       */
+      template <int dim, int spacedim, int structdim>
+      static types::global_dof_index
+      get_dof_index(const DoFHandler<dim, spacedim>              &dof_handler,
+                    const unsigned int                            obj_level,
+                    const unsigned int                            obj_index,
+                    const types::fe_index                         fe_index,
+                    const unsigned int                            local_index,
+                    const std::integral_constant<int, structdim> &dd)
+      {
+        types::global_dof_index global_index;
+        process_dof_index(dof_handler,
+                          obj_level,
+                          obj_index,
+                          fe_index,
+                          local_index,
+                          dd,
+                          global_index,
+                          [](const auto &ptr, auto &value) { value = ptr; });
+        return global_index;
+      }
+
+
+      template <int dim, int spacedim>
+      static types::global_dof_index &
+      mg_vertex_dof_index(DoFHandler<dim, spacedim> &dof_handler,
+                          const int                  level,
+                          const unsigned int         vertex_index,
+                          const unsigned int         i)
+      {
+        Assert(dof_handler.hp_capability_enabled == false,
+               ExcMessage(
+                 "DoFHandler in hp-mode does not implement multilevel DoFs."));
+
+        return dof_handler.mg_vertex_dofs[vertex_index].access_index(
+          level, i, dof_handler.get_fe().n_dofs_per_vertex());
+      }
+
+
+
+      /**
+       * Return the number of different finite elements that are active on a
+       * given object such as a vertex, line, or cell.
+       *
+       * The template argument `structdim` indicates the
+       * dimensionality of the object on which we seek to know the DoF
+       * index. For example, if `structdim==0`, then we are looking to
+       * get a DoF index on a vertex of the indicated cell.
+       */
+      template <int dim, int spacedim, int structdim>
+      static unsigned int
+      n_active_fe_indices(const DoFHandler<dim, spacedim> &dof_handler,
+                          const unsigned int               obj_level,
+                          const unsigned int               obj_index,
+                          const std::integral_constant<int, structdim> &)
+      {
+        (void)obj_level;
+
+        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
+
+        // 1) no hp used -> fe_index == 0
+        if (dof_handler.hp_capability_enabled == false)
+          return 1;
+
+        // 2) cell and hp is used -> there is only one fe_index
+        if (structdim == dim)
+          return 1;
+
+        // 3) general entity and hp is used
+        AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
+        AssertIndexRange(obj_index + 1,
+                         dof_handler.hp_object_fe_ptr[structdim].size());
+
+        return dof_handler.hp_object_fe_ptr[structdim][obj_index + 1] -
+               dof_handler.hp_object_fe_ptr[structdim][obj_index];
+      }
+
+
+
+      /**
+       * Return the FE index of the n-th finite element active on a given
+       * object such as a vertex, line, or cell.
+       *
+       * The template argument `structdim` indicates the
+       * dimensionality of the object on which we seek to know the DoF
+       * index. For example, if `structdim==0`, then we are looking to
+       * get a DoF index on a vertex of the indicated cell.
+       */
+      template <int dim, int spacedim, int structdim>
+      static types::fe_index
+      nth_active_fe_index(const DoFHandler<dim, spacedim> &dof_handler,
+                          const unsigned int               obj_level,
+                          const unsigned int               obj_index,
+                          const unsigned int               local_index,
+                          const std::integral_constant<int, structdim> &)
+      {
+        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
+
+        // for cells only one active FE index available
+        Assert(((structdim == dim) &&
+                (local_index != DoFHandler<dim, spacedim>::default_fe_index)) ==
+                 false,
+               ExcNotImplemented());
+
+        // 1) no hp used -> fe_index == 0
+        if (dof_handler.hp_capability_enabled == false)
+          return DoFHandler<dim, spacedim>::default_fe_index;
+
+        // 2) cell and hp is used -> there is only one fe_index
+        if (structdim == dim)
+          return dof_handler.hp_cell_active_fe_indices[obj_level][obj_index];
+
+        // 3) general entity and hp is used
+        AssertIndexRange(structdim, dof_handler.hp_object_fe_indices.size());
+        AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
+        AssertIndexRange(obj_index,
+                         dof_handler.hp_object_fe_ptr[structdim].size());
+        AssertIndexRange(dof_handler.hp_object_fe_ptr[structdim][obj_index] +
+                           local_index,
+                         dof_handler.hp_object_fe_indices[structdim].size());
+
+        return dof_handler.hp_object_fe_indices
+          [structdim]
+          [dof_handler.hp_object_fe_ptr[structdim][obj_index] + local_index];
+      }
+
+
+
+      /**
+       * Returns all active FE indices on a given object such as a
+       * vertex, line, or cell.
+       *
+       * The size of the returned set equals the number of finite elements that
+       * are active on this vertex.
+       *
+       * The template argument `structdim` indicates the
+       * dimensionality of the object on which we seek to know the DoF
+       * index. For example, if `structdim==0`, then we are looking to
+       * get a DoF index on a vertex of the indicated cell.
+       */
+      template <int dim, int spacedim, int structdim>
+      static std::set<types::fe_index>
+      get_active_fe_indices(const DoFHandler<dim, spacedim> &dof_handler,
+                            const unsigned int               obj_level,
+                            const unsigned int               obj_index,
+                            const std::integral_constant<int, structdim> &t)
+      {
+        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
+
+        // 1) no hp used -> fe_index == 0
+        if (dof_handler.hp_capability_enabled == false)
+          return {DoFHandler<dim, spacedim>::default_fe_index};
+
+        // 2) cell and hp is used -> there is only one fe_index
+        if (structdim == dim)
+          return {dof_handler.hp_cell_active_fe_indices[obj_level][obj_index]};
+
+        // 3) general entity and hp is used
+        std::set<types::fe_index> active_fe_indices;
+        for (unsigned int i = 0;
+             i < n_active_fe_indices(dof_handler, obj_level, obj_index, t);
+             ++i)
+          active_fe_indices.insert(
+            nth_active_fe_index(dof_handler, obj_level, obj_index, i, t));
+        return active_fe_indices;
+      }
+
+
+
+      template <int dim, int spacedim, int structdim>
+      static bool
+      fe_index_is_active(const DoFHandler<dim, spacedim> &dof_handler,
+                         const unsigned int               obj_level,
+                         const unsigned int               obj_index,
+                         const types::fe_index            fe_index,
+                         const std::integral_constant<int, structdim> &)
+      {
+        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
+
+        // 1) no hp used -> fe_index == 0
+        if (dof_handler.hp_capability_enabled == false)
+          return (fe_index == DoFHandler<dim, spacedim>::default_fe_index);
+
+        // 2) cell and hp is used -> there is only one fe_index
+        if (structdim == dim)
+          return dof_handler.hp_cell_active_fe_indices[obj_level][obj_index] ==
+                 fe_index;
+
+        // 3) general entity and hp is used
+        return std::find(
+                 dof_handler.hp_object_fe_indices[structdim].begin() +
+                   dof_handler.hp_object_fe_ptr[structdim][obj_index],
+                 dof_handler.hp_object_fe_indices[structdim].begin() +
+                   dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
+                 fe_index) !=
+               (dof_handler.hp_object_fe_indices[structdim].begin() +
+                dof_handler.hp_object_fe_ptr[structdim][obj_index + 1]);
+      }
+
+
+
+      template <typename InputVector, typename ForwardIterator>
+      static void
+      extract_subvector_to(const InputVector             &values,
+                           const types::global_dof_index *cache,
+                           const types::global_dof_index *cache_end,
+                           ForwardIterator                local_values_begin)
+      {
+        values.extract_subvector_to(cache, cache_end, local_values_begin);
+      }
+
+
+
+#ifdef DEAL_II_WITH_TRILINOS
+      static std::vector<unsigned int>
+      sort_indices(const types::global_dof_index *v_begin,
+                   const types::global_dof_index *v_end)
+      {
+        // initialize original index locations
+        std::vector<unsigned int> idx(v_end - v_begin);
+        std::iota(idx.begin(), idx.end(), 0u);
+
+        // sort indices based on comparing values in v
+        std::sort(idx.begin(),
+                  idx.end(),
+                  [&v_begin](unsigned int i1, unsigned int i2) {
+                    return *(v_begin + i1) < *(v_begin + i2);
+                  });
+
+        return idx;
+      }
+
+
+
+#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
+      template <typename ForwardIterator, typename Number, typename MemorySpace>
+      static void
+      extract_subvector_to(
+        const LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace>
+                                      &values,
+        const types::global_dof_index *cache_begin,
+        const types::global_dof_index *cache_end,
+        ForwardIterator                local_values_begin)
+      {
+        std::vector<unsigned int> sorted_indices_pos =
+          sort_indices(cache_begin, cache_end);
+        const unsigned int cache_size = cache_end - cache_begin;
+        std::vector<types::global_dof_index> cache_indices(cache_size);
+        for (unsigned int i = 0; i < cache_size; ++i)
+          cache_indices[i] = *(cache_begin + sorted_indices_pos[i]);
+
+        IndexSet index_set(cache_indices.back() + 1);
+        index_set.add_indices(cache_indices.begin(), cache_indices.end());
+        index_set.compress();
+        LinearAlgebra::ReadWriteVector<Number> read_write_vector(index_set);
+        read_write_vector.import_elements(values, VectorOperation::insert);
+
+        // Copy the elements from read_write_vector and reorder them.
+        for (unsigned int i = 0; i < cache_size; ++i, ++local_values_begin)
+          *local_values_begin = read_write_vector[sorted_indices_pos[i]];
+      }
+#  endif
+
+
+
+      template <typename ForwardIterator>
+      static void
+      extract_subvector_to(const LinearAlgebra::EpetraWrappers::Vector &values,
+                           const types::global_dof_index *cache_begin,
+                           const types::global_dof_index *cache_end,
+                           ForwardIterator                local_values_begin)
+      {
+        std::vector<unsigned int> sorted_indices_pos =
+          sort_indices(cache_begin, cache_end);
+        const unsigned int cache_size = cache_end - cache_begin;
+        std::vector<types::global_dof_index> cache_indices(cache_size);
+        for (unsigned int i = 0; i < cache_size; ++i)
+          cache_indices[i] = *(cache_begin + sorted_indices_pos[i]);
+
+        IndexSet index_set(cache_indices.back() + 1);
+        index_set.add_indices(cache_indices.begin(), cache_indices.end());
+        index_set.compress();
+        LinearAlgebra::ReadWriteVector<double> read_write_vector(index_set);
+        read_write_vector.import_elements(values, VectorOperation::insert);
+
+        // Copy the elements from read_write_vector and reorder them.
+        for (unsigned int i = 0; i < cache_size; ++i, ++local_values_begin)
+          *local_values_begin = read_write_vector[sorted_indices_pos[i]];
+      }
+#endif
+
+      /**
+       * Loop over all degrees of freedom of the object described by the
+       * provided @p accessor and @p fe_index and count them.
+       */
+      template <int dim, int spacedim, bool level_dof_access, int structdim>
+      static unsigned int
+      n_dof_indices(
+        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
+                             &accessor,
+        const types::fe_index fe_index_,
+        const bool            count_level_dofs)
+      {
+        // note: we cannot rely on the template parameter level_dof_access here,
+        // since the function get_mg_dof_indices()/set_mg_dof_indices() can be
+        // called even if level_dof_access==false.
+        if (count_level_dofs)
+          {
+            const auto &fe = accessor.get_fe(fe_index_);
+
+            const unsigned int                                   //
+              dofs_per_vertex = fe.n_dofs_per_vertex(),          //
+              dofs_per_line   = fe.n_dofs_per_line(),            //
+              dofs_per_quad   = fe.n_dofs_per_quad(0 /*dummy*/), //
+              dofs_per_hex    = fe.n_dofs_per_hex();             //
+
+            unsigned int index = 0;
+
+            // 1) VERTEX dofs
+            index += dofs_per_vertex * accessor.n_vertices();
+
+            // 2) LINE dofs
+            if (structdim == 2 || structdim == 3)
+              index += dofs_per_line * accessor.n_lines();
+
+            // 3) FACE dofs
+            if (structdim == 3)
+              index += dofs_per_quad * accessor.n_faces();
+
+            // 4) INNER dofs
+            const unsigned int interior_dofs =
+              structdim == 1 ? dofs_per_line :
+                               (structdim == 2 ? dofs_per_quad : dofs_per_hex);
+
+            index += interior_dofs;
+
+            return index;
+          }
+        else
+          {
+            const auto fe_index =
+              internal::DoFAccessorImplementation::get_fe_index_or_default(
+                accessor, fe_index_);
+
+            unsigned int index = 0;
+
+            // 1) VERTEX dofs
+            for (const auto vertex : accessor.vertex_indices())
+              index += process_object_range(accessor.get_dof_handler(),
+                                            0,
+                                            accessor.vertex_index(vertex),
+                                            fe_index,
+                                            std::integral_constant<int, 0>())
+                         .second;
+
+            // 2) LINE dofs
+            if (structdim == 2 || structdim == 3)
+              for (const auto line : accessor.line_indices())
+                index +=
+                  process_object_range(*accessor.line(line), fe_index).second;
+
+            // 3) FACE dofs
+            if (structdim == 3)
+              for (const auto face : accessor.face_indices())
+                index +=
+                  process_object_range(*accessor.quad(face), fe_index).second;
+
+            // 4) INNER dofs
+            index += process_object_range(accessor, fe_index).second;
+
+            return index;
+          }
+      }
+
+
+
+      // The next few internal helper functions are needed to support various
+      // DoFIndicesType kinds, e.g. actual vectors of DoFIndices or empty
+      // types that we use when we only want to work on the internally stored
+      // DoFs and never extract any number.
+      template <typename ArrayType>
+      static unsigned int
+      get_array_length(const ArrayType &array)
+      {
+        return array.size();
+      }
+
+      static unsigned int
+      get_array_length(const std::tuple<> &)
+      {
+        return 0;
+      }
+
+      template <typename ArrayType>
+      static types::global_dof_index *
+      get_array_ptr(const ArrayType &array)
+      {
+        return const_cast<types::global_dof_index *>(array.data());
+      }
+
+      static types::global_dof_index *
+      get_array_ptr(const std::tuple<> &)
+      {
+        return nullptr;
+      }
+
+
+
+      /**
+       * Loop over all degrees of freedom of the object described by the
+       * provided @p accessor and @p fe_index and perform the static functions
+       * provided by DoFOperation (set/get) on these.
+       */
+      template <int  dim,
+                int  spacedim,
+                bool level_dof_access,
+                int  structdim,
+                typename DoFIndicesType,
+                typename DoFOperation,
+                typename DoFProcessor>
+      static void
+      process_dof_indices(
+        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
+                             &accessor,
+        const DoFIndicesType &const_dof_indices,
+        const types::fe_index fe_index_,
+        const DoFOperation   &dof_operation,
+        const DoFProcessor   &dof_processor,
+        const bool            count_level_dofs)
+      {
+        const types::fe_index fe_index =
+          internal::DoFAccessorImplementation::get_fe_index_or_default(
+            accessor, fe_index_);
+
+        // we cannot rely on the template parameter level_dof_access here, since
+        // the function get_mg_dof_indices()/set_mg_dof_indices() can be called
+        // even if level_dof_access==false.
+        (void)count_level_dofs;
+
+        const auto &fe = accessor.get_fe(fe_index);
+
+        // we want to pass in rvalue 'std::tuple<>' types as `DoFIndicesType`,
+        // but we need non-const references for std::vector<> types, so get in
+        // a const reference here and immediately cast the constness away -
+        // note that any use of the dereferenced invalid type will result in a
+        // segfault
+        types::global_dof_index *dof_indices_ptr =
+          get_array_ptr(const_dof_indices);
+        types::global_dof_index *end_dof_indices =
+          dof_indices_ptr + get_array_length(const_dof_indices);
+
+        // 1) VERTEX dofs, only step into the functions if we actually have
+        // DoFs on them
+        if (fe.n_dofs_per_vertex() > 0)
+          for (const auto vertex : accessor.vertex_indices())
+            dof_operation.process_vertex_dofs(*accessor.dof_handler,
+                                              accessor.vertex_index(vertex),
+                                              fe_index,
+                                              dof_indices_ptr,
+                                              dof_processor);
+
+        // 2) copy dof numbers from the LINE. for lines with the wrong
+        // orientation (which might occur in 3d), we have already made sure that
+        // we're ok by picking the correct vertices (this happens automatically
+        // in the vertex() function). however, if the line is in wrong
+        // orientation, we look at it in flipped orientation and we will have to
+        // adjust the shape function indices that we see to correspond to the
+        // correct (face/cell-local) ordering.
+        if ((structdim == 2 || structdim == 3) && fe.n_dofs_per_line() > 0)
+          {
+            const auto line_indices = internal::TriaAccessorImplementation::
+              Implementation::get_line_indices_of_cell(accessor);
+            const auto line_orientations =
+              internal::TriaAccessorImplementation::Implementation::
+                get_line_orientations_of_cell(accessor);
+
+            for (const auto line : accessor.line_indices())
+              {
+                const bool line_orientation = line_orientations[line];
+                if (line_orientation)
+                  dof_operation.process_dofs(
+                    accessor.get_dof_handler(),
+                    0,
+                    line_indices[line],
+                    fe_index,
+                    [](const auto d) { return d; },
+                    std::integral_constant<int, 1>(),
+                    dof_indices_ptr,
+                    dof_processor);
+                else
+                  dof_operation.process_dofs(
+                    accessor.get_dof_handler(),
+                    0,
+                    line_indices[line],
+                    fe_index,
+                    [&fe, line_orientation](const auto d) {
+                      return fe.adjust_line_dof_index_for_line_orientation(
+                        d, line_orientation);
+                    },
+                    std::integral_constant<int, 1>(),
+                    dof_indices_ptr,
+                    dof_processor);
+              }
+          }
+
+        // 3) copy dof numbers from the FACE. for faces with the wrong
+        // orientation, we have already made sure that we're ok by picking the
+        // correct lines and vertices (this happens automatically in the line()
+        // and vertex() functions). however, if the face is in wrong
+        // orientation, we look at it in flipped orientation and we will have to
+        // adjust the shape function indices that we see to correspond to the
+        // correct (cell-local) ordering. The same applies, if the face_rotation
+        // or face_orientation is non-standard
+        if (structdim == 3 && fe.max_dofs_per_quad() > 0)
+          for (const auto face_no : accessor.face_indices())
+            {
+              const auto combined_orientation =
+                accessor.combined_face_orientation(face_no);
+              const unsigned int quad_index = accessor.quad_index(face_no);
+              if (combined_orientation ==
+                  ReferenceCell::default_combined_face_orientation())
+                dof_operation.process_dofs(
+                  accessor.get_dof_handler(),
+                  0,
+                  quad_index,
+                  fe_index,
+                  [](const auto d) { return d; },
+                  std::integral_constant<int, 2>(),
+                  dof_indices_ptr,
+                  dof_processor);
+              else
+                dof_operation.process_dofs(
+                  accessor.get_dof_handler(),
+                  0,
+                  quad_index,
+                  fe_index,
+                  [&](const auto d) {
+                    return fe.adjust_quad_dof_index_for_face_orientation(
+                      d, face_no, combined_orientation);
+                  },
+                  std::integral_constant<int, 2>(),
+                  dof_indices_ptr,
+                  dof_processor);
+            }
+
+        // 4) INNER dofs - here we need to make sure that the shortcut to not
+        // run the function does not miss the faces of wedge and pyramid
+        // elements where n_dofs_per_object might not return the largest
+        // possible value
+        if (((dim == 3 && structdim == 2) ?
+               fe.max_dofs_per_quad() :
+               fe.template n_dofs_per_object<structdim>()) > 0)
+          dof_operation.process_dofs(
+            accessor.get_dof_handler(),
+            accessor.level(),
+            accessor.index(),
+            fe_index,
+            [&](const auto d) { return d; },
+            std::integral_constant<int, structdim>(),
+            dof_indices_ptr,
+            dof_processor);
+
+        if (dof_indices_ptr != nullptr)
+          {
+            AssertDimension(n_dof_indices(accessor, fe_index, count_level_dofs),
+                            dof_indices_ptr - get_array_ptr(const_dof_indices));
+          }
+
+        // PM: This is a part that should not be reached since it indicates that
+        // an object (and/or its subobjects) is not active. However,
+        // unfortunately this function is called by
+        // DoFTools::set_periodicity_constraints() indirectly by
+        // get_dof_indices() also for artificial faces to determine if a face
+        // is artificial.
+        types::global_dof_index invalid_index = numbers::invalid_dof_index;
+        for (; dof_indices_ptr < end_dof_indices; ++dof_indices_ptr)
+          dof_processor(invalid_index, dof_indices_ptr);
+      }
+
+
+
+      /**
+       * An internal struct encapsulating the task of getting (vertex)
+       * DoF indices.
+       */
+      template <int dim, int spacedim>
+      struct DoFIndexProcessor
+      {
+        /**
+         * Return vertex DoF indices.
+         */
+        template <typename DoFProcessor>
+        DEAL_II_ALWAYS_INLINE void
+        process_vertex_dofs(DoFHandler<dim, spacedim> &dof_handler,
+                            const unsigned int         vertex_index,
+                            const types::fe_index      fe_index,
+                            types::global_dof_index  *&dof_indices_ptr,
+                            const DoFProcessor        &dof_processor) const
+        {
+          process_object(
+            dof_handler,
+            0,
+            vertex_index,
+            fe_index,
+            [](const auto d) {
+              DEAL_II_ASSERT_UNREACHABLE();
+              return d;
+            },
+            std::integral_constant<int, 0>(),
+            dof_indices_ptr,
+            dof_processor);
+        }
+
+        /**
+         * Return DoF indices for lines, quads, and inner degrees of freedom.
+         */
+        template <int structdim, typename DoFMapping, typename DoFProcessor>
+        DEAL_II_ALWAYS_INLINE void
+        process_dofs(const DoFHandler<dim, spacedim> &dof_handler,
+                     const unsigned int               obj_level,
+                     const unsigned int               obj_index,
+                     const types::fe_index            fe_index,
+                     const DoFMapping                &mapping,
+                     const std::integral_constant<int, structdim>,
+                     types::global_dof_index *&dof_indices_ptr,
+                     const DoFProcessor       &dof_processor) const
+        {
+          process_object(
+            dof_handler,
+            obj_level,
+            obj_index,
+            fe_index,
+            mapping,
+            std::integral_constant<int, std::min(structdim, dim)>(),
+            dof_indices_ptr,
+            dof_processor);
+        }
+      };
+
+
+
+      /**
+       * An internal struct encapsulating the task of getting level (vertex)
+       * DoF indices.
+       */
+      template <int dim, int spacedim>
+      struct MGDoFIndexProcessor
+      {
+        /**
+         * Constructor.
+         */
+        MGDoFIndexProcessor(const unsigned int level)
+          : level(level)
+        {}
+
+        /**
+         * Return vertex DoF indices.
+         */
+        template <typename DoFProcessor>
+        DEAL_II_ALWAYS_INLINE void
+        process_vertex_dofs(DoFHandler<dim, spacedim> &dof_handler,
+                            const unsigned int         vertex_index,
+                            const types::fe_index,
+                            types::global_dof_index *&dof_indices_ptr,
+                            const DoFProcessor       &dof_processor) const
+        {
+          const unsigned int n_indices =
+            dof_handler.get_fe(0).template n_dofs_per_object<0>();
+          types::global_dof_index *stored_indices =
+            &dof_handler.mg_vertex_dofs[vertex_index].access_index(level,
+                                                                   0,
+                                                                   n_indices);
+          for (unsigned int d = 0; d < n_indices; ++d, ++dof_indices_ptr)
+            dof_processor(stored_indices[d], dof_indices_ptr);
+        }
+
+        /**
+         * Return DoF indices for lines, quads, and inner degrees of freedom.
+         */
+        template <int structdim, typename DoFMapping, typename DoFProcessor>
+        DEAL_II_ALWAYS_INLINE void
+        process_dofs(const DoFHandler<dim, spacedim> &dof_handler,
+                     const unsigned int,
+                     const unsigned int    obj_index,
+                     const types::fe_index fe_index,
+                     const DoFMapping     &mapping,
+                     const std::integral_constant<int, structdim>,
+                     types::global_dof_index *&dof_indices_ptr,
+                     const DoFProcessor       &dof_processor) const
+        {
+          const unsigned int n_indices =
+            dof_handler.get_fe(0).template n_dofs_per_object<structdim>();
+          types::global_dof_index *stored_indices = &get_mg_dof_index(
+            dof_handler,
+            dof_handler.mg_levels[level],
+            dof_handler.mg_faces,
+            obj_index,
+            fe_index,
+            0,
+            std::integral_constant<int, std::min(structdim, dim)>());
+          for (unsigned int d = 0; d < n_indices; ++d, ++dof_indices_ptr)
+            dof_processor(stored_indices[structdim < dim ? mapping(d) : d],
+                          dof_indices_ptr);
+        }
+
+      private:
+        const unsigned int level;
+      };
+
+
+
+      template <int dim, int spacedim, bool level_dof_access, int structdim>
+      static void
+      get_dof_indices(
+        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
+                                             &accessor,
+        std::vector<types::global_dof_index> &dof_indices,
+        const types::fe_index                 fe_index)
+      {
+        process_dof_indices(
+          accessor,
+          dof_indices,
+          fe_index,
+          DoFIndexProcessor<dim, spacedim>(),
+          [](auto stored_index, auto dof_ptr) { *dof_ptr = stored_index; },
+          false);
+      }
+
+
+
+      template <int dim, int spacedim, bool level_dof_access, int structdim>
+      static void
+      set_dof_indices(
+        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
+                                                   &accessor,
+        const std::vector<types::global_dof_index> &dof_indices,
+        const types::fe_index                       fe_index)
+      {
+        // Note: this function is as general as `get_dof_indices()`. This
+        // assert is placed here since it is currently only used by the
+        // function DoFCellAccessor::set_dof_indices(), which is called by
+        // internal::DoFHandlerImplementation::Policy::Implementation::distribute_dofs().
+        // In the case of new use cases, this assert can be removed.
+        Assert(
+          dim == structdim,
+          ExcMessage(
+            "This function is intended to be used for DoFCellAccessor, i.e., "
+            "dimension == structdim."));
+
+        process_dof_indices(
+          accessor,
+          dof_indices,
+          fe_index,
+          DoFIndexProcessor<dim, spacedim>(),
+          [](auto &stored_index, auto dof_ptr) { stored_index = *dof_ptr; },
+          false);
+      }
+
+
+
+      template <int dim, int spacedim, bool level_dof_access, int structdim>
+      static void
+      get_mg_dof_indices(
+        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
+                                             &accessor,
+        const int                             level,
+        std::vector<types::global_dof_index> &dof_indices,
+        const types::fe_index                 fe_index)
+      {
+        Assert((fe_index == DoFHandler<dim, spacedim>::default_fe_index),
+               ExcMessage("MG DoF indices cannot be queried in hp case"));
+        process_dof_indices(
+          accessor,
+          dof_indices,
+          fe_index,
+          MGDoFIndexProcessor<dim, spacedim>(level),
+          [](auto stored_index, auto dof_ptr) { *dof_ptr = stored_index; },
+          true);
+      }
+
+
+
+      template <int dim, int spacedim, bool level_dof_access, int structdim>
+      static void
+      set_mg_dof_indices(
+        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
+                                                   &accessor,
+        const int                                   level,
+        const std::vector<types::global_dof_index> &dof_indices,
+        const types::fe_index                       fe_index)
+      {
+        Assert((fe_index == DoFHandler<dim, spacedim>::default_fe_index),
+               ExcMessage("MG DoF indices cannot be queried in hp case"));
+
+        // Note: this function is as general as `get_mg_dof_indices()`. This
+        // assert is placed here since it is currently only used by the
+        // function DoFCellAccessor::set_mg_dof_indices(), which is called by
+        // internal::DoFHandlerImplementation::Policy::Implementation::distribute_mg_dofs().
+        // In the case of new use cases, this assert can be removed.
+        Assert(dim == structdim,
+               ExcMessage("This function is intended to be used for "
+                          "DoFCellAccessor, i.e., dimension == structdim."));
+
+        process_dof_indices(
+          accessor,
+          dof_indices,
+          fe_index,
+          MGDoFIndexProcessor<dim, spacedim>(level),
+          [](auto &stored_index, auto dof_ptr) { stored_index = *dof_ptr; },
+          true);
+      }
+
+
+
+      template <int dim, int spacedim>
+      static types::global_dof_index &
+      get_mg_dof_index(
+        const DoFHandler<dim, spacedim> &dof_handler,
+        const std::unique_ptr<internal::DoFHandlerImplementation::DoFLevel<dim>>
+          &mg_level,
+        const std::unique_ptr<internal::DoFHandlerImplementation::DoFFaces<dim>>
+          &,
+        const unsigned int    obj_index,
+        const types::fe_index fe_index,
+        const unsigned int    local_index,
+        const std::integral_constant<int, dim>)
+      {
+        Assert(dof_handler.hp_capability_enabled == false,
+               (typename DoFHandler<dim, spacedim>::ExcNotImplementedWithHP()));
+
+        return mg_level->dof_object.access_dof_index(
+          static_cast<const DoFHandler<dim, spacedim> &>(dof_handler),
+          obj_index,
+          fe_index,
+          local_index);
+      }
+
+
+
+      template <int dim, int spacedim, std::enable_if_t<(dim > 1), int> = 0>
+      static types::global_dof_index &
+      get_mg_dof_index(
+        const DoFHandler<dim, spacedim> &dof_handler,
+        const std::unique_ptr<internal::DoFHandlerImplementation::DoFLevel<dim>>
+          &,
+        const std::unique_ptr<internal::DoFHandlerImplementation::DoFFaces<dim>>
+                             &mg_faces,
+        const unsigned int    obj_index,
+        const types::fe_index fe_index,
+        const unsigned int    local_index,
+        const std::integral_constant<int, 1>)
+      {
+        return mg_faces->lines.access_dof_index(
+          static_cast<const DoFHandler<dim, spacedim> &>(dof_handler),
+          obj_index,
+          fe_index,
+          local_index);
+      }
+
+
+
+      template <int spacedim>
+      static types::global_dof_index &
+      get_mg_dof_index(
+        const DoFHandler<3, spacedim> &dof_handler,
+        const std::unique_ptr<internal::DoFHandlerImplementation::DoFLevel<3>>
+          &,
+        const std::unique_ptr<internal::DoFHandlerImplementation::DoFFaces<3>>
+                             &mg_faces,
+        const unsigned int    obj_index,
+        const types::fe_index fe_index,
+        const unsigned int    local_index,
+        const std::integral_constant<int, 2>)
+      {
+        Assert(dof_handler.hp_capability_enabled == false,
+               (typename DoFHandler<3, spacedim>::ExcNotImplementedWithHP()));
+        return mg_faces->quads.access_dof_index(
+          static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
+          obj_index,
+          fe_index,
+          local_index);
+      }
+    };
+
+
+
+    template <int dim, int spacedim, bool level_dof_access>
+    void
+    get_cell_dof_indices(
+      const dealii::DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
+      Implementation::dof_index_vector_type &dof_indices,
+      const unsigned int                     fe_index);
+  } // namespace DoFAccessorImplementation
+} // namespace internal
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline types::global_dof_index
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::dof_index(
+  const unsigned int    i,
+  const types::fe_index fe_index_) const
+{
+  const auto fe_index =
+    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
+                                                                 fe_index_);
+
+  // access the respective DoF
+  return dealii::internal::DoFAccessorImplementation::Implementation::
+    get_dof_index(*this->dof_handler,
+                  this->level(),
+                  this->index(),
+                  fe_index,
+                  i,
+                  std::integral_constant<int, structdim>());
+}
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline types::global_dof_index
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::mg_dof_index(
+  const int          level,
+  const unsigned int i) const
+{
+  return internal::DoFAccessorImplementation::Implementation::get_mg_dof_index(
+    *this->dof_handler,
+    this->dof_handler->mg_levels[level],
+    this->dof_handler->mg_faces,
+    this->index(),
+    0,
+    i,
+    std::integral_constant<int, structdim>());
+}
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline void
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::set_dof_index(
+  const unsigned int            i,
+  const types::global_dof_index index,
+  const types::fe_index         fe_index_) const
+{
+  const auto fe_index =
+    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
+                                                                 fe_index_);
+
+  // access the respective DoF
+  dealii::internal::DoFAccessorImplementation::Implementation::set_dof_index(
+    *this->dof_handler,
+    this->level(),
+    this->index(),
+    fe_index,
+    i,
+    std::integral_constant<int, structdim>(),
+    index);
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline void
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::set_mg_dof_index(
+  const int                     level,
+  const unsigned int            i,
+  const types::global_dof_index index) const
+{
+  internal::DoFAccessorImplementation::Implementation::get_mg_dof_index(
+    *this->dof_handler,
+    this->dof_handler->mg_levels[level],
+    this->dof_handler->mg_faces,
+    this->index(),
+    0,
+    i,
+    std::integral_constant<int, structdim>()) = index;
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline unsigned int
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::n_active_fe_indices()
+  const
+{
+  // access the respective DoF
+  return dealii::internal::DoFAccessorImplementation::Implementation::
+    n_active_fe_indices(*this->dof_handler,
+                        this->level(),
+                        this->index(),
+                        std::integral_constant<int, structdim>());
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline types::fe_index
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::nth_active_fe_index(
+  const unsigned int n) const
+{
+  // access the respective DoF
+  return dealii::internal::DoFAccessorImplementation::Implementation::
+    nth_active_fe_index(*this->dof_handler,
+                        this->level(),
+                        this->index(),
+                        n,
+                        std::integral_constant<int, structdim>());
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline std::set<types::fe_index>
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::get_active_fe_indices()
+  const
+{
+  std::set<types::fe_index> active_fe_indices;
+  for (unsigned int i = 0; i < n_active_fe_indices(); ++i)
+    active_fe_indices.insert(nth_active_fe_index(i));
+  return active_fe_indices;
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline bool
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::fe_index_is_active(
+  const types::fe_index fe_index) const
+{
+  // access the respective DoF
+  return dealii::internal::DoFAccessorImplementation::Implementation::
+    fe_index_is_active(*this->dof_handler,
+                       this->level(),
+                       this->index(),
+                       fe_index,
+                       std::integral_constant<int, structdim>());
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline types::global_dof_index
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::vertex_dof_index(
+  const unsigned int    vertex,
+  const unsigned int    i,
+  const types::fe_index fe_index_) const
+{
+  const types::fe_index fe_index =
+    (((this->dof_handler->hp_capability_enabled == false) &&
+      (fe_index_ == numbers::invalid_fe_index)) ?
+       // No hp enabled, and the argument is at its default value -> we
+       // can translate to the default active fe index
+       DoFHandler<dim, spacedim>::default_fe_index :
+       // Otherwise: If anything other than the default is provided by
+       // the caller, then we should take just that. As an exception, if
+       // we are on a cell (rather than a face/edge/vertex), then we know
+       // that there is only one active fe index on this cell and we can
+       // use that:
+       ((dim == structdim) && (fe_index_ == numbers::invalid_fe_index) ?
+          this->nth_active_fe_index(0) :
+          fe_index_));
+
+  return dealii::internal::DoFAccessorImplementation::Implementation::
+    get_dof_index(*this->dof_handler,
+                  0,
+                  this->vertex_index(vertex),
+                  fe_index,
+                  i,
+                  std::integral_constant<int, 0>());
+}
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline types::global_dof_index
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::mg_vertex_dof_index(
+  const int             level,
+  const unsigned int    vertex,
+  const unsigned int    i,
+  const types::fe_index fe_index_) const
+{
+  const auto fe_index =
+    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
+                                                                 fe_index_);
+  (void)fe_index;
+  Assert(this->dof_handler != nullptr, ExcInvalidObject());
+  Assert(this->dof_handler->mg_vertex_dofs.size() > 0,
+         ExcMessage("Multigrid DoF indices can only be accessed after "
+                    "DoFHandler::distribute_mg_dofs() has been called!"));
+  AssertIndexRange(vertex, this->n_vertices());
+  AssertIndexRange(i, this->dof_handler->get_fe(fe_index).n_dofs_per_vertex());
+
+  Assert(dof_handler->hp_capability_enabled == false,
+         ExcMessage(
+           "DoFHandler in hp-mode does not implement multilevel DoFs."));
+
+  return this->dof_handler->mg_vertex_dofs[this->vertex_index(vertex)]
+    .access_index(level, i, this->dof_handler->get_fe().n_dofs_per_vertex());
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline void
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::
+  set_mg_vertex_dof_index(const int                     level,
+                          const unsigned int            vertex,
+                          const unsigned int            i,
+                          const types::global_dof_index index,
+                          const types::fe_index         fe_index_) const
+{
+  const auto fe_index =
+    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
+                                                                 fe_index_);
+  (void)fe_index;
+  Assert(this->dof_handler != nullptr, ExcInvalidObject());
+  AssertIndexRange(vertex, this->n_vertices());
+  AssertIndexRange(i, this->dof_handler->get_fe(fe_index).n_dofs_per_vertex());
+
+  Assert(dof_handler->hp_capability_enabled == false,
+         ExcMessage(
+           "DoFHandler in hp-mode does not implement multilevel DoFs."));
+
+  this->dof_handler->mg_vertex_dofs[this->vertex_index(vertex)].access_index(
+    level, i, this->dof_handler->get_fe().n_dofs_per_vertex()) = index;
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline const FiniteElement<dim, spacedim> &
+DoFAccessor<structdim, dim, spacedim, level_dof_access>::get_fe(
+  const types::fe_index fe_index) const
+{
+  Assert(fe_index_is_active(fe_index) == true,
+         ExcMessage("This function can only be called for active FE indices"));
+
+  return this->dof_handler->get_fe(fe_index);
+}
+
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline typename dealii::internal::DoFHandlerImplementation::
+  Iterators<dim, spacedim, level_dof_access>::line_iterator
+  DoFAccessor<structdim, dim, spacedim, level_dof_access>::line(
+    const unsigned int i) const
+{
+  // if we are asking for a particular line and this object refers to
+  // a line, then the only valid index is i==0 and we should return
+  // *this
+  if (structdim == 1)
+    {
+      Assert(i == 0,
+             ExcMessage("You can only ask for line zero if the "
+                        "current object is a line itself."));
+      return typename dealii::internal::DoFHandlerImplementation::
+        Iterators<dim, spacedim, level_dof_access>::cell_iterator(
+          &this->get_triangulation(),
+          this->level(),
+          this->index(),
+          &this->get_dof_handler());
+    }
+
+  // otherwise we need to be in structdim>=2
+  Assert(structdim > 1, ExcImpossibleInDim(structdim));
+  Assert(dim > 1, ExcImpossibleInDim(dim));
+
+  // checking of 'i' happens in line_index(i)
+  return typename dealii::internal::DoFHandlerImplementation::
+    Iterators<dim, spacedim, level_dof_access>::line_iterator(
+      this->tria,
+      0, // only sub-objects are allowed, which have no level
+      this->line_index(i),
+      this->dof_handler);
+}
+
+
+template <int structdim, int dim, int spacedim, bool level_dof_access>
+inline typename dealii::internal::DoFHandlerImplementation::
+  Iterators<dim, spacedim, level_dof_access>::quad_iterator
+  DoFAccessor<structdim, dim, spacedim, level_dof_access>::quad(
+    const unsigned int i) const
+{
+  // if we are asking for a
+  // particular quad and this object
+  // refers to a quad, then the only
+  // valid index is i==0 and we
+  // should return *this
+  if (structdim == 2)
+    {
+      Assert(i == 0,
+             ExcMessage("You can only ask for quad zero if the "
+                        "current object is a quad itself."));
+      return typename dealii::internal::DoFHandlerImplementation::
+        Iterators<dim, spacedim>::cell_iterator(&this->get_triangulation(),
+                                                this->level(),
+                                                this->index(),
+                                                &this->get_dof_handler());
+    }
+
+  // otherwise we need to be in structdim>=3
+  Assert(structdim > 2, ExcImpossibleInDim(structdim));
+  Assert(dim > 2, ExcImpossibleInDim(dim));
+
+  // checking of 'i' happens in quad_index(i)
+  return typename dealii::internal::DoFHandlerImplementation::
+    Iterators<dim, spacedim, level_dof_access>::quad_iterator(
+      this->tria,
+      0, // only sub-objects are allowed, which have no level
+      this->quad_index(i),
+      this->dof_handler);
+}
+
+
+/*----------------- Functions: DoFAccessor<0,1,spacedim> --------------------*/
+
+
+template <int spacedim, bool level_dof_access>
+inline DoFAccessor<0, 1, spacedim, level_dof_access>::DoFAccessor()
+{
+  Assert(false, ExcInvalidObject());
+}
+
+
 
-// include more templates
-#include <deal.II/dofs/dof_accessor.templates.h>
+template <int spacedim, bool level_dof_access>
+inline DoFAccessor<0, 1, spacedim, level_dof_access>::DoFAccessor(
+  const Triangulation<1, spacedim>                       *tria,
+  const typename TriaAccessor<0, 1, spacedim>::VertexKind vertex_kind,
+  const unsigned int                                      vertex_index,
+  const DoFHandler<1, spacedim>                          *dof_handler)
+  : BaseClass(tria, vertex_kind, vertex_index)
+  , dof_handler(const_cast<DoFHandler<1, spacedim> *>(dof_handler))
+{}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline DoFAccessor<0, 1, spacedim, level_dof_access>::DoFAccessor(
+  const Triangulation<1, spacedim> *tria,
+  const int                         level,
+  const int                         index,
+  const DoFHandler<1, spacedim>    *dof_handler)
+  // This is the constructor signature for "ordinary" (non-vertex)
+  // accessors and we shouldn't be calling it altogether. But it is also
+  // the constructor that the default-constructor of TriaRawIterator
+  // calls when default-constructing an iterator object. If so, this
+  // happens with level==-2 and index==-2, and this is the only case we
+  // would like to support. We do this by just forwarding to the
+  // other constructor of this class, and then asserting the condition
+  // on level and index.
+  : DoFAccessor<0, 1, spacedim, level_dof_access>(
+      tria,
+      TriaAccessor<0, 1, spacedim>::interior_vertex,
+      0U,
+      dof_handler)
+{
+  (void)level;
+  (void)index;
+  Assert((tria == nullptr) && (level == -2) && (index == -2) &&
+           (dof_handler == nullptr),
+         ExcMessage(
+           "This constructor can not be called for face iterators in 1d, "
+           "except to default-construct iterator objects."));
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+template <int structdim2, int dim2, int spacedim2>
+DoFAccessor<0, 1, spacedim, level_dof_access>::DoFAccessor(
+  const InvalidAccessor<structdim2, dim2, spacedim2> &)
+{
+  Assert(false, ExcInvalidObject());
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
+inline DoFAccessor<0, 1, spacedim, level_dof_access>::DoFAccessor(
+  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &)
+{
+  Assert(false, ExcInvalidObject());
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline void
+DoFAccessor<0, 1, spacedim, level_dof_access>::set_dof_handler(
+  DoFHandler<1, spacedim> *dh)
+{
+  Assert(dh != nullptr, ExcInvalidObject());
+  this->dof_handler = dh;
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline void
+DoFAccessor<0, 1, spacedim, level_dof_access>::set_dof_index(
+  const unsigned int /*i*/,
+  const types::global_dof_index /*index*/,
+  const types::fe_index /*fe_index*/) const
+{
+  DEAL_II_NOT_IMPLEMENTED();
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline const DoFHandler<1, spacedim> &
+DoFAccessor<0, 1, spacedim, level_dof_access>::get_dof_handler() const
+{
+  return *this->dof_handler;
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline void
+DoFAccessor<0, 1, spacedim, level_dof_access>::get_dof_indices(
+  std::vector<types::global_dof_index> &dof_indices,
+  const types::fe_index                 fe_index_) const
+{
+  const auto fe_index =
+    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
+                                                                 fe_index_);
+
+  for (unsigned int i = 0; i < dof_indices.size(); ++i)
+    dof_indices[i] = dealii::internal::DoFAccessorImplementation::
+      Implementation::get_dof_index(*dof_handler,
+                                    0,
+                                    this->global_vertex_index,
+                                    fe_index,
+                                    i,
+                                    std::integral_constant<int, 0>());
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline void
+DoFAccessor<0, 1, spacedim, level_dof_access>::get_mg_dof_indices(
+  const int                             level,
+  std::vector<types::global_dof_index> &dof_indices,
+  const types::fe_index                 fe_index_) const
+{
+  const auto fe_index =
+    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
+                                                                 fe_index_);
+  (void)fe_index;
+  AssertDimension(fe_index, (DoFHandler<1, spacedim>::default_fe_index));
+
+  for (unsigned int i = 0; i < dof_indices.size(); ++i)
+    dof_indices[i] =
+      dealii::internal::DoFAccessorImplementation::Implementation::
+        mg_vertex_dof_index(*dof_handler, level, this->global_vertex_index, i);
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline types::global_dof_index
+DoFAccessor<0, 1, spacedim, level_dof_access>::vertex_dof_index(
+  const unsigned int    vertex,
+  const unsigned int    i,
+  const types::fe_index fe_index_) const
+{
+  const auto fe_index =
+    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
+                                                                 fe_index_);
+
+  (void)vertex;
+  AssertIndexRange(vertex, 1);
+  return dealii::internal::DoFAccessorImplementation::Implementation::
+    get_dof_index(*dof_handler,
+                  0,
+                  this->global_vertex_index,
+                  fe_index,
+                  i,
+                  std::integral_constant<int, 0>());
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline types::global_dof_index
+DoFAccessor<0, 1, spacedim, level_dof_access>::dof_index(
+  const unsigned int    i,
+  const types::fe_index fe_index_) const
+{
+  const auto fe_index =
+    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
+                                                                 fe_index_);
+
+  return dealii::internal::DoFAccessorImplementation::Implementation::
+    get_dof_index(*this->dof_handler,
+                  0,
+                  this->vertex_index(0),
+                  fe_index,
+                  i,
+                  std::integral_constant<int, 0>());
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline unsigned int
+DoFAccessor<0, 1, spacedim, level_dof_access>::n_active_fe_indices() const
+{
+  return 1;
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline types::fe_index
+DoFAccessor<0, 1, spacedim, level_dof_access>::nth_active_fe_index(
+  const unsigned int /*n*/) const
+{
+  return 0;
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline bool
+DoFAccessor<0, 1, spacedim, level_dof_access>::fe_index_is_active(
+  const types::fe_index /*fe_index*/) const
+{
+  DEAL_II_NOT_IMPLEMENTED();
+  return false;
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline const FiniteElement<1, spacedim> &
+DoFAccessor<0, 1, spacedim, level_dof_access>::get_fe(
+  const types::fe_index fe_index) const
+{
+  Assert(this->dof_handler != nullptr, ExcInvalidObject());
+  return dof_handler->get_fe(fe_index);
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline void
+DoFAccessor<0, 1, spacedim, level_dof_access>::copy_from(
+  const TriaAccessorBase<0, 1, spacedim> &da)
+{
+  Assert(this->dof_handler != nullptr, ExcInvalidObject());
+  BaseClass::copy_from(da);
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+template <bool level_dof_access2>
+inline void
+DoFAccessor<0, 1, spacedim, level_dof_access>::copy_from(
+  const DoFAccessor<0, 1, spacedim, level_dof_access2> &a)
+{
+  BaseClass::copy_from(a);
+  set_dof_handler(a.dof_handler);
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline TriaIterator<DoFAccessor<0, 1, spacedim, level_dof_access>>
+DoFAccessor<0, 1, spacedim, level_dof_access>::child(
+  const unsigned int /*i*/) const
+{
+  return TriaIterator<DoFAccessor<0, 1, spacedim, level_dof_access>>();
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline typename dealii::internal::DoFHandlerImplementation::
+  Iterators<1, spacedim, level_dof_access>::line_iterator
+  DoFAccessor<0, 1, spacedim, level_dof_access>::line(
+    const unsigned int /*c*/) const
+{
+  DEAL_II_NOT_IMPLEMENTED();
+  return typename dealii::internal::DoFHandlerImplementation::
+    Iterators<1, spacedim, level_dof_access>::line_iterator();
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+inline typename dealii::internal::DoFHandlerImplementation::
+  Iterators<1, spacedim, level_dof_access>::quad_iterator
+  DoFAccessor<0, 1, spacedim, level_dof_access>::quad(
+    const unsigned int /*c*/) const
+{
+  DEAL_II_NOT_IMPLEMENTED();
+  return typename dealii::internal::DoFHandlerImplementation::
+    Iterators<1, spacedim, level_dof_access>::quad_iterator();
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
+inline bool
+DoFAccessor<0, 1, spacedim, level_dof_access>::operator==(
+  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &a) const
+{
+  Assert(structdim2 == 0, ExcCantCompareIterators());
+  Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
+  return (BaseClass::operator==(a));
+}
+
+
+
+template <int spacedim, bool level_dof_access>
+template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
+inline bool
+DoFAccessor<0, 1, spacedim, level_dof_access>::operator!=(
+  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &a) const
+{
+  Assert(structdim2 == 0, ExcCantCompareIterators());
+  Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
+  return (BaseClass::operator!=(a));
+}
+
+
+
+/*------------------------- Functions: DoFCellAccessor -----------------------*/
+
+
+namespace internal
+{
+  namespace DoFCellAccessorImplementation
+  {
+    /**
+     * A class with the same purpose as the similarly named class of the
+     * Triangulation class. See there for more information.
+     */
+    struct Implementation
+    {
+      /**
+       * Do what the active_fe_index function in the parent class is supposed to
+       * do.
+       */
+      template <int dim, int spacedim, bool level_dof_access>
+      static types::fe_index
+      active_fe_index(
+        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor)
+      {
+        if (accessor.dof_handler->hp_capability_enabled == false)
+          return DoFHandler<dim, spacedim>::default_fe_index;
+
+        Assert(accessor.dof_handler != nullptr,
+               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
+        Assert(static_cast<unsigned int>(accessor.level()) <
+                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
+               ExcMessage("DoFHandler not initialized"));
+
+        return accessor.dof_handler
+          ->hp_cell_active_fe_indices[accessor.level()][accessor.present_index];
+      }
+
+
+
+      /**
+       * Do what the set_active_fe_index function in the parent class is
+       * supposed to do.
+       */
+      template <int dim, int spacedim, bool level_dof_access>
+      static void
+      set_active_fe_index(
+        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
+        const types::fe_index                                   i)
+      {
+        if (accessor.dof_handler->hp_capability_enabled == false)
+          {
+            AssertDimension(i, (DoFHandler<dim, spacedim>::default_fe_index));
+            return;
+          }
+
+        Assert(accessor.dof_handler != nullptr,
+               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
+        Assert(static_cast<unsigned int>(accessor.level()) <
+                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
+               ExcMessage("DoFHandler not initialized"));
+        Assert(i != numbers::invalid_fe_index,
+               ExcMessage("Invalid finite element index."));
+
+        accessor.dof_handler
+          ->hp_cell_active_fe_indices[accessor.level()]
+                                     [accessor.present_index] = i;
+      }
+
+
+
+      /**
+       * Do what the future_fe_index function in the parent class is supposed to
+       * do.
+       */
+      template <int dim, int spacedim, bool level_dof_access>
+      static types::fe_index
+      future_fe_index(
+        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor)
+      {
+        if (accessor.dof_handler->hp_capability_enabled == false)
+          return DoFHandler<dim, spacedim>::default_fe_index;
+
+        Assert(accessor.dof_handler != nullptr,
+               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
+        Assert(static_cast<unsigned int>(accessor.level()) <
+                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
+               ExcMessage("DoFHandler not initialized"));
+
+        if (future_fe_index_set(accessor))
+          return accessor.dof_handler
+            ->hp_cell_future_fe_indices[accessor.level()]
+                                       [accessor.present_index];
+        else
+          return accessor.dof_handler
+            ->hp_cell_active_fe_indices[accessor.level()]
+                                       [accessor.present_index];
+      }
+
+
+      /**
+       * Do what the set_future_fe_index function in the parent class is
+       * supposed to do.
+       */
+      template <int dim, int spacedim, bool level_dof_access>
+      static void
+      set_future_fe_index(
+        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
+        const types::fe_index                                   i)
+      {
+        if (accessor.dof_handler->hp_capability_enabled == false)
+          {
+            AssertDimension(i, (DoFHandler<dim, spacedim>::default_fe_index));
+            return;
+          }
+
+        Assert(accessor.dof_handler != nullptr,
+               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
+        Assert(static_cast<unsigned int>(accessor.level()) <
+                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
+               ExcMessage("DoFHandler not initialized"));
+        Assert(i != numbers::invalid_fe_index,
+               ExcMessage("Invalid finite element index."));
+
+        accessor.dof_handler
+          ->hp_cell_future_fe_indices[accessor.level()]
+                                     [accessor.present_index] = i;
+      }
+
+
+
+      /**
+       * Do what the future_fe_index_set function in the parent class is
+       * supposed to do.
+       */
+      template <int dim, int spacedim, bool level_dof_access>
+      static bool
+      future_fe_index_set(
+        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor)
+      {
+        if (accessor.dof_handler->hp_capability_enabled == false)
+          return false;
+
+        Assert(accessor.dof_handler != nullptr,
+               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
+        Assert(static_cast<unsigned int>(accessor.level()) <
+                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
+               ExcMessage("DoFHandler not initialized"));
+
+        return accessor.dof_handler
+                 ->hp_cell_future_fe_indices[accessor.level()]
+                                            [accessor.present_index] !=
+               numbers::invalid_fe_index;
+      }
+
+
+
+      /**
+       * Do what the clear_fe_index function in the parent class is supposed to
+       * do.
+       */
+      template <int dim, int spacedim, bool level_dof_access>
+      static void
+      clear_future_fe_index(
+        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor)
+      {
+        if (accessor.dof_handler->hp_capability_enabled == false)
+          return;
+
+        Assert(accessor.dof_handler != nullptr,
+               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
+        Assert(static_cast<unsigned int>(accessor.level()) <
+                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
+               ExcMessage("DoFHandler not initialized"));
+
+        accessor.dof_handler
+          ->hp_cell_future_fe_indices[accessor.level()]
+                                     [accessor.present_index] =
+          numbers::invalid_fe_index;
+      }
+    };
+  } // namespace DoFCellAccessorImplementation
+} // namespace internal
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  DoFCellAccessor(const Triangulation<dimension_, space_dimension_> *tria,
+                  const int                                          level,
+                  const int                                          index,
+                  const AccessorData                                *local_data)
+  : DoFAccessor<dimension_, dimension_, space_dimension_, level_dof_access>(
+      tria,
+      level,
+      index,
+      local_data)
+{}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <int structdim2, int dim2, int spacedim2>
+inline DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  DoFCellAccessor(const InvalidAccessor<structdim2, dim2, spacedim2> &)
+{
+  Assert(false, typename BaseClass::ExcInvalidObject());
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
+inline DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  DoFCellAccessor(
+    const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &other)
+  : BaseClass(other)
+{}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline TriaIterator<
+  DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::neighbor(
+  const unsigned int i) const
+{
+  TriaIterator<DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
+    q(this->tria,
+      this->neighbor_level(i),
+      this->neighbor_index(i),
+      this->dof_handler);
+
+#ifdef DEBUG
+  if (q.state() != IteratorState::past_the_end)
+    Assert(q->used(), ExcInternalError());
+#endif
+  return q;
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline TriaIterator<
+  DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::child(
+  const unsigned int i) const
+{
+  TriaIterator<DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
+    q(this->tria, this->level() + 1, this->child_index(i), this->dof_handler);
+
+#ifdef DEBUG
+  if (q.state() != IteratorState::past_the_end)
+    Assert(q->used(), ExcInternalError());
+#endif
+  return q;
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline boost::container::small_vector<
+  TriaIterator<DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>,
+  GeometryInfo<dimension_>::max_children_per_cell>
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  child_iterators() const
+{
+  boost::container::small_vector<
+    TriaIterator<
+      DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>,
+    GeometryInfo<dimension_>::max_children_per_cell>
+    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 dimension_, int space_dimension_, bool level_dof_access>
+inline TriaIterator<
+  DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::parent() const
+{
+  TriaIterator<DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
+    q(this->tria, this->level() - 1, this->parent_index(), this->dof_handler);
+
+  return q;
+}
+
+
+
+namespace internal
+{
+  namespace DoFCellAccessorImplementation
+  {
+    template <int dim, int spacedim, bool level_dof_access>
+    inline TriaIterator<
+      dealii::DoFAccessor<dim - 1, dim, spacedim, level_dof_access>>
+    get_face(
+      const dealii::DoFCellAccessor<dim, spacedim, level_dof_access> &cell,
+      const unsigned int                                              i,
+      const std::integral_constant<int, 1>)
+    {
+      dealii::DoFAccessor<0, dim, spacedim, level_dof_access> a(
+        &cell.get_triangulation(),
+        ((i == 0) && cell.at_boundary(0) ?
+           dealii::TriaAccessor<0, 1, spacedim>::left_vertex :
+           ((i == 1) && cell.at_boundary(1) ?
+              dealii::TriaAccessor<0, 1, spacedim>::right_vertex :
+              dealii::TriaAccessor<0, 1, spacedim>::interior_vertex)),
+        cell.vertex_index(i),
+        &cell.get_dof_handler());
+      return dealii::TriaIterator<
+        dealii::DoFAccessor<0, dim, spacedim, level_dof_access>>(a);
+    }
+
+
+    template <int dim, int spacedim, bool level_dof_access>
+    inline TriaIterator<
+      dealii::DoFAccessor<dim - 1, dim, spacedim, level_dof_access>>
+    get_face(
+      const dealii::DoFCellAccessor<dim, spacedim, level_dof_access> &cell,
+      const unsigned int                                              i,
+      const std::integral_constant<int, 2>)
+    {
+      return cell.line(i);
+    }
+
+
+    template <int dim, int spacedim, bool level_dof_access>
+    inline TriaIterator<
+      dealii::DoFAccessor<dim - 1, dim, spacedim, level_dof_access>>
+    get_face(
+      const dealii::DoFCellAccessor<dim, spacedim, level_dof_access> &cell,
+      const unsigned int                                              i,
+      const std::integral_constant<int, 3>)
+    {
+      return cell.quad(i);
+    }
+  } // namespace DoFCellAccessorImplementation
+} // namespace internal
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline typename DoFCellAccessor<dimension_,
+                                space_dimension_,
+                                level_dof_access>::face_iterator
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::face(
+  const unsigned int i) const
+{
+  AssertIndexRange(i, this->n_faces());
+
+  return dealii::internal::DoFCellAccessorImplementation::get_face(
+    *this, i, std::integral_constant<int, dimension_>());
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline boost::container::small_vector<
+  typename DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+    face_iterator,
+  GeometryInfo<dimension_>::faces_per_cell>
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  face_iterators() const
+{
+  boost::container::small_vector<
+    typename DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+      face_iterator,
+    GeometryInfo<dimension_>::faces_per_cell>
+    face_iterators(this->n_faces());
+
+  for (const unsigned int i : this->face_indices())
+    face_iterators[i] =
+      dealii::internal::DoFCellAccessorImplementation::get_face(
+        *this, i, std::integral_constant<int, dimension_>());
+
+  return face_iterators;
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  get_active_or_mg_dof_indices(
+    std::vector<types::global_dof_index> &dof_indices) const
+{
+  if (level_dof_access)
+    get_mg_dof_indices(dof_indices);
+  else
+    get_dof_indices(dof_indices);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <class InputVector, typename number>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::get_dof_values(
+  const InputVector &values,
+  Vector<number>    &local_values) const
+{
+  get_dof_values(values, local_values.begin(), local_values.end());
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <typename Number, typename ForwardIterator>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::get_dof_values(
+  const ReadVector<Number> &values,
+  ForwardIterator           local_values_begin,
+  ForwardIterator           local_values_end) const
+{
+  (void)local_values_end;
+  Assert(this->is_artificial() == false,
+         ExcMessage("Can't ask for DoF indices on artificial cells."));
+  Assert(this->is_active(), ExcMessage("Cell must be active."));
+  Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
+
+  Assert(static_cast<unsigned int>(local_values_end - local_values_begin) ==
+           this->get_fe().n_dofs_per_cell(),
+         typename DoFCellAccessor::ExcVectorDoesNotMatch());
+  Assert(values.size() == this->get_dof_handler().n_dofs(),
+         typename DoFCellAccessor::ExcVectorDoesNotMatch());
+
+  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
+    dof_indices(this->get_fe().n_dofs_per_cell());
+  internal::DoFAccessorImplementation::get_cell_dof_indices(
+    *this, dof_indices, this->active_fe_index());
+
+  boost::container::small_vector<Number, 27> values_temp(local_values_end -
+                                                         local_values_begin);
+  auto view = make_array_view(values_temp.begin(), values_temp.end());
+  values.extract_subvector_to(make_array_view(dof_indices.begin(),
+                                              dof_indices.end()),
+                              view);
+  using view_type = std::remove_reference_t<decltype(*local_values_begin)>;
+  ArrayView<view_type> values_view2(&*local_values_begin,
+                                    local_values_end - local_values_begin);
+  std::copy(values_temp.begin(), values_temp.end(), values_view2.begin());
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <class InputVector, typename ForwardIterator>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::get_dof_values(
+  const AffineConstraints<typename InputVector::value_type> &constraints,
+  const InputVector                                         &values,
+  ForwardIterator                                            local_values_begin,
+  ForwardIterator local_values_end) const
+{
+  Assert(this->is_artificial() == false,
+         ExcMessage("Can't ask for DoF indices on artificial cells."));
+  Assert(this->is_active(), ExcMessage("Cell must be active."));
+
+  Assert(static_cast<unsigned int>(local_values_end - local_values_begin) ==
+           this->get_fe().n_dofs_per_cell(),
+         typename DoFCellAccessor::ExcVectorDoesNotMatch());
+  Assert(values.size() == this->get_dof_handler().n_dofs(),
+         typename DoFCellAccessor::ExcVectorDoesNotMatch());
+
+
+  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
+    dof_indices(this->get_fe().n_dofs_per_cell());
+  internal::DoFAccessorImplementation::get_cell_dof_indices(
+    *this, dof_indices, this->active_fe_index());
+
+  constraints.get_dof_values(values,
+                             dof_indices.data(),
+                             local_values_begin,
+                             local_values_end);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <class OutputVector, typename number>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::set_dof_values(
+  const Vector<number> &local_values,
+  OutputVector         &values) const
+{
+  Assert(this->is_artificial() == false,
+         ExcMessage("Can't ask for DoF indices on artificial cells."));
+  Assert(this->is_active(), ExcMessage("Cell must be active."));
+
+  Assert(static_cast<unsigned int>(local_values.size()) ==
+           this->get_fe().n_dofs_per_cell(),
+         typename DoFCellAccessor::ExcVectorDoesNotMatch());
+  Assert(values.size() == this->get_dof_handler().n_dofs(),
+         typename DoFCellAccessor::ExcVectorDoesNotMatch());
+
+
+  Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
+  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
+    dof_indices(this->get_fe().n_dofs_per_cell());
+  internal::DoFAccessorImplementation::get_cell_dof_indices(
+    *this, dof_indices, this->active_fe_index());
+
+  for (unsigned int i = 0; i < this->get_fe().n_dofs_per_cell(); ++i)
+    internal::ElementAccess<OutputVector>::set(local_values(i),
+                                               dof_indices[i],
+                                               values);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline const FiniteElement<dimension_, space_dimension_> &
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::get_fe() const
+{
+  Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           this->is_active(),
+         ExcMessage(
+           "For DoFHandler objects in hp-mode, finite elements are only "
+           "associated with active cells. Consequently, you can not ask "
+           "for the active finite element on cells with children."));
+
+  const auto &fe = this->dof_handler->get_fe(active_fe_index());
+
+  Assert(this->reference_cell() == fe.reference_cell(),
+         internal::ExcNonMatchingReferenceCellTypes(this->reference_cell(),
+                                                    fe.reference_cell()));
+
+  return fe;
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline types::fe_index
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  active_fe_index() const
+{
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           this->is_active(),
+         ExcMessage(
+           "You can not ask for the active FE index on a cell that has "
+           "children because no degrees of freedom are assigned "
+           "to this cell and, consequently, no finite element "
+           "is associated with it."));
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           (this->is_locally_owned() || this->is_ghost()),
+         ExcMessage("You can only query active FE index information on cells "
+                    "that are either locally owned or (after distributing "
+                    "degrees of freedom) are ghost cells."));
+
+  return dealii::internal::DoFCellAccessorImplementation::Implementation::
+    active_fe_index(*this);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  set_active_fe_index(const types::fe_index i) const
+{
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           this->is_active(),
+         ExcMessage("You can not set the active FE index on a cell that has "
+                    "children because no degrees of freedom will be assigned "
+                    "to this cell."));
+
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           this->is_locally_owned(),
+         ExcMessage("You can only set active FE index information on cells "
+                    "that are locally owned. On ghost cells, this information "
+                    "will automatically be propagated from the owning process "
+                    "of that cell, and there is no information at all on "
+                    "artificial cells."));
+
+  dealii::internal::DoFCellAccessorImplementation::Implementation::
+    set_active_fe_index(*this, i);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline const FiniteElement<dimension_, space_dimension_> &
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::get_future_fe()
+  const
+{
+  Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           this->is_active(),
+         ExcMessage(
+           "For DoFHandler objects in hp-mode, finite elements are only "
+           "associated with active cells. Consequently, you can not ask "
+           "for the future finite element on cells with children."));
+
+  return this->dof_handler->get_fe(future_fe_index());
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline types::fe_index
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  future_fe_index() const
+{
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           (this->has_children() == false),
+         ExcMessage(
+           "You can not ask for the future FE index on a cell that has "
+           "children because no degrees of freedom are assigned "
+           "to this cell and, consequently, no finite element "
+           "is associated with it."));
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           (this->is_locally_owned()),
+         ExcMessage("You can only query future FE index information on cells "
+                    "that are locally owned."));
+
+  return dealii::internal::DoFCellAccessorImplementation::Implementation::
+    future_fe_index(*this);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  set_future_fe_index(const types::fe_index i) const
+{
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           (this->has_children() == false),
+         ExcMessage("You can not set the future FE index on a cell that has "
+                    "children because no degrees of freedom will be assigned "
+                    "to this cell."));
+
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           this->is_locally_owned(),
+         ExcMessage("You can only set future FE index information on cells "
+                    "that are locally owned."));
+
+  dealii::internal::DoFCellAccessorImplementation::Implementation::
+    set_future_fe_index(*this, i);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline bool
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  future_fe_index_set() const
+{
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           (this->has_children() == false),
+         ExcMessage(
+           "You can not ask for the future FE index on a cell that has "
+           "children because no degrees of freedom are assigned "
+           "to this cell and, consequently, no finite element "
+           "is associated with it."));
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           (this->is_locally_owned()),
+         ExcMessage("You can only query future FE index information on cells "
+                    "that are locally owned."));
+
+  return dealii::internal::DoFCellAccessorImplementation::Implementation::
+    future_fe_index_set(*this);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  clear_future_fe_index() const
+{
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           (this->has_children() == false),
+         ExcMessage(
+           "You can not ask for the future FE index on a cell that has "
+           "children because no degrees of freedom are assigned "
+           "to this cell and, consequently, no finite element "
+           "is associated with it."));
+  Assert((this->dof_handler->hp_capability_enabled == false) ||
+           (this->is_locally_owned()),
+         ExcMessage("You can only query future FE index information on cells "
+                    "that are locally owned."));
+
+  dealii::internal::DoFCellAccessorImplementation::Implementation::
+    clear_future_fe_index(*this);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <typename number, typename OutputVector>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  distribute_local_to_global(const Vector<number> &local_source,
+                             OutputVector         &global_destination) const
+{
+  this->distribute_local_to_global(local_source.begin(),
+                                   local_source.end(),
+                                   global_destination);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <typename ForwardIterator, typename OutputVector>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  distribute_local_to_global(ForwardIterator local_source_begin,
+                             ForwardIterator local_source_end,
+                             OutputVector   &global_destination) const
+{
+  Assert(this->dof_handler != nullptr,
+         (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
+  Assert(static_cast<unsigned int>(local_source_end - local_source_begin) ==
+           this->get_fe().n_dofs_per_cell(),
+         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
+  Assert(this->dof_handler->n_dofs() == global_destination.size(),
+         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
+
+  Assert(!this->has_children(), ExcMessage("Cell must be active"));
+
+  Assert(
+    internal::ArrayViewHelper::is_contiguous(local_source_begin,
+                                             local_source_end),
+    ExcMessage(
+      "This function can not be called with iterator types that do not point to contiguous memory."));
+
+  const unsigned int n_dofs = local_source_end - local_source_begin;
+
+  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
+    dof_indices(n_dofs);
+  internal::DoFAccessorImplementation::get_cell_dof_indices(
+    *this, dof_indices, this->active_fe_index());
+
+  // distribute cell vector
+  global_destination.add(n_dofs, dof_indices.data(), &(*local_source_begin));
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <typename ForwardIterator, typename OutputVector>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  distribute_local_to_global(
+    const AffineConstraints<typename OutputVector::value_type> &constraints,
+    ForwardIterator local_source_begin,
+    ForwardIterator local_source_end,
+    OutputVector   &global_destination) const
+{
+  Assert(this->dof_handler != nullptr,
+         (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
+  Assert(local_source_end - local_source_begin ==
+           this->get_fe().n_dofs_per_cell(),
+         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
+  Assert(this->dof_handler->n_dofs() == global_destination.size(),
+         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
+
+  Assert(!this->has_children(), ExcMessage("Cell must be active."));
+
+  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
+    dof_indices(this->get_fe().n_dofs_per_cell());
+  internal::DoFAccessorImplementation::get_cell_dof_indices(
+    *this, dof_indices, this->active_fe_index());
+
+  // distribute cell vector
+  constraints.distribute_local_to_global(local_source_begin,
+                                         local_source_end,
+                                         dof_indices.data(),
+                                         global_destination);
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <typename number, typename OutputMatrix>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  distribute_local_to_global(const FullMatrix<number> &local_source,
+                             OutputMatrix             &global_destination) const
+{
+  Assert(this->dof_handler != nullptr,
+         (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
+  Assert(local_source.m() == this->get_fe().n_dofs_per_cell(),
+         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
+  Assert(local_source.n() == this->get_fe().n_dofs_per_cell(),
+         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
+  Assert(this->dof_handler->n_dofs() == global_destination.m(),
+         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
+  Assert(this->dof_handler->n_dofs() == global_destination.n(),
+         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
+
+  Assert(!this->has_children(), ExcMessage("Cell must be active."));
+
+  const unsigned int n_dofs = local_source.m();
+
+  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
+    dof_indices(n_dofs);
+  internal::DoFAccessorImplementation::get_cell_dof_indices(
+    *this, dof_indices, this->active_fe_index());
+
+  // distribute cell matrix
+  for (unsigned int i = 0; i < n_dofs; ++i)
+    global_destination.add(dof_indices[i],
+                           n_dofs,
+                           dof_indices.data(),
+                           &local_source(i, 0));
+}
+
+
+
+template <int dimension_, int space_dimension_, bool level_dof_access>
+template <typename number, typename OutputMatrix, typename OutputVector>
+inline void
+DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
+  distribute_local_to_global(const FullMatrix<number> &local_matrix,
+                             const Vector<number>     &local_vector,
+                             OutputMatrix             &global_matrix,
+                             OutputVector             &global_vector) const
+{
+  Assert(this->dof_handler != nullptr,
+         (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
+  Assert(local_matrix.m() == this->get_fe().n_dofs_per_cell(),
+         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
+  Assert(local_matrix.n() == this->get_fe().n_dofs_per_cell(),
+         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
+  Assert(this->dof_handler->n_dofs() == global_matrix.m(),
+         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
+  Assert(this->dof_handler->n_dofs() == global_matrix.n(),
+         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
+  Assert(local_vector.size() == this->get_fe().n_dofs_per_cell(),
+         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
+  Assert(this->dof_handler->n_dofs() == global_vector.size(),
+         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
+
+  Assert(!this->has_children(), ExcMessage("Cell must be active."));
+
+  const unsigned int n_dofs = this->get_fe().n_dofs_per_cell();
+  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
+    dof_indices(n_dofs);
+  internal::DoFAccessorImplementation::get_cell_dof_indices(
+    *this, dof_indices, this->active_fe_index());
+
+  // distribute cell matrices
+  for (unsigned int i = 0; i < n_dofs; ++i)
+    {
+      global_matrix.add(dof_indices[i],
+                        n_dofs,
+                        dof_indices.data(),
+                        &local_matrix(i, 0));
+      global_vector(dof_indices[i]) += local_vector(i);
+    }
+}
+
+
+
+DEAL_II_NAMESPACE_CLOSE
 
 
 #endif
index 851b3e1bd43c7584987dd6a9ae485516bdd3710f..c544aa5f3d7500e2b392f58f7debd0d8d223d667 100644 (file)
 #ifndef dealii_dof_accessor_templates_h
 #define dealii_dof_accessor_templates_h
 
-
-#include <deal.II/base/config.h>
-
-#include <deal.II/base/types.h>
-
-#include <deal.II/dofs/dof_accessor.h>
-#include <deal.II/dofs/dof_faces.h>
-#include <deal.II/dofs/dof_levels.h>
-
-#include <deal.II/grid/tria_iterator.h>
-#include <deal.II/grid/tria_iterator.templates.h>
-
-#include <deal.II/lac/affine_constraints.h>
-#include <deal.II/lac/read_write_vector.h>
-
-#include <boost/container/small_vector.hpp>
-
-#include <limits>
-#include <type_traits>
-#include <vector>
-
-DEAL_II_NAMESPACE_OPEN
-
-
-/*------------------------- Functions: DoFAccessor ---------------------------*/
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline DoFAccessor<structdim, dim, spacedim, level_dof_access>::DoFAccessor()
-{
-  Assert(false, ExcInvalidObject());
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline DoFAccessor<structdim, dim, spacedim, level_dof_access>::DoFAccessor(
-  const Triangulation<dim, spacedim> *tria,
-  const int                           level,
-  const int                           index,
-  const DoFHandler<dim, spacedim>    *dof_handler)
-  : dealii::internal::DoFAccessorImplementation::
-      Inheritance<structdim, dim, spacedim>::BaseClass(tria, level, index)
-  , dof_handler(const_cast<DoFHandler<dim, spacedim> *>(dof_handler))
-{
-  Assert(
-    tria == nullptr || &dof_handler->get_triangulation() == tria,
-    ExcMessage(
-      "You can't create a DoF accessor in which the DoFHandler object "
-      "uses a different triangulation than the one you pass as argument."));
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-template <int structdim2, int dim2, int spacedim2>
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::DoFAccessor(
-  const InvalidAccessor<structdim2, dim2, spacedim2> &)
-{
-  Assert(false, ExcInvalidObject());
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
-inline DoFAccessor<structdim, dim, spacedim, level_dof_access>::DoFAccessor(
-  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &other)
-  : BaseClass(other)
-  , dof_handler(nullptr)
-{
-  Assert(false,
-         ExcMessage(
-           "You are trying to assign iterators that are incompatible. "
-           "The reason for incompatibility is that they refer to objects of "
-           "different dimensionality (e.g., assigning a line iterator "
-           "to a quad iterator)."));
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-template <bool level_dof_access2>
-inline DoFAccessor<structdim, dim, spacedim, level_dof_access>::DoFAccessor(
-  const DoFAccessor<structdim, dim, spacedim, level_dof_access2> &other)
-  : BaseClass(other)
-  , dof_handler(const_cast<DoFHandler<dim, spacedim> *>(other.dof_handler))
-{}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline void
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::set_dof_handler(
-  DoFHandler<dim, spacedim> *dh)
-{
-  Assert(dh != nullptr, ExcInvalidObject());
-  this->dof_handler = dh;
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline const DoFHandler<dim, spacedim> &
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::get_dof_handler() const
-{
-  Assert(this->dof_handler != nullptr, ExcInvalidObject());
-  return *this->dof_handler;
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline void
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::copy_from(
-  const TriaAccessorBase<structdim, dim, spacedim> &da)
-{
-  Assert(this->dof_handler != nullptr, ExcInvalidObject());
-  BaseClass::copy_from(da);
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-template <bool level_dof_access2>
-inline void
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::copy_from(
-  const DoFAccessor<structdim, dim, spacedim, level_dof_access2> &a)
-{
-  BaseClass::copy_from(a);
-  this->dof_handler = a.dof_handler;
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
-inline bool
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::operator==(
-  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &a) const
-{
-  Assert(structdim == structdim2, ExcCantCompareIterators());
-  Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
-  return (BaseClass::operator==(a));
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
-inline bool
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::operator!=(
-  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &a) const
-{
-  Assert(structdim == structdim2, ExcCantCompareIterators());
-  Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
-  return (BaseClass::operator!=(a));
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline TriaIterator<DoFAccessor<structdim, dim, spacedim, level_dof_access>>
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::child(
-  const unsigned int i) const
-{
-  Assert(static_cast<unsigned int>(this->present_level) <
-           this->dof_handler->object_dof_indices.size(),
-         ExcMessage("DoFHandler not initialized"));
-
-  TriaIterator<TriaAccessor<structdim, dim, spacedim>> t =
-    TriaAccessor<structdim, dim, spacedim>::child(i);
-
-  TriaIterator<DoFAccessor<structdim, dim, spacedim, level_dof_access>> q(
-    *t, this->dof_handler);
-  return q;
-}
-
-
-namespace internal
-{
-  namespace DoFAccessorImplementation
-  {
-    /**
-     * Convert an FE index that might contain the right value but also
-     * numbers::invalid_fe_index to a right value if needed/possible.
-     */
-    template <int structdim, int dim, int spacedim, bool level_dof_access>
-    types::fe_index
-    get_fe_index_or_default(
-      const DoFAccessor<structdim, dim, spacedim, level_dof_access> &cell,
-      const types::fe_index                                          fe_index)
-    {
-      if (cell.get_dof_handler().has_hp_capabilities() == false)
-        {
-          // No hp enabled, and the argument is at its default value -> we
-          // can translate to the default active fe index
-          Assert(
-            (fe_index == numbers::invalid_fe_index) ||
-              (fe_index == DoFHandler<dim, spacedim>::default_fe_index),
-            ExcMessage(
-              "It is not possible to specify a FE index if no hp support is used!"));
-
-          return DoFHandler<dim, spacedim>::default_fe_index;
-        }
-      else
-        {
-          // Otherwise: If anything other than the default is provided by
-          // the caller, then we should take just that. As an exception, if
-          // we are on a cell (rather than a face/edge/vertex), then we know
-          // that there is only one active fe index on this cell and we can
-          // use that:
-          if ((dim == structdim) && (fe_index == numbers::invalid_fe_index))
-            {
-              AssertDimension(cell.n_active_fe_indices(), 1);
-
-              return cell.nth_active_fe_index(0);
-            }
-
-          Assert((fe_index != numbers::invalid_fe_index),
-                 ExcMessage(
-                   "You need to specify a FE index if hp support is used!"));
-
-          return fe_index;
-        }
-    }
-
-    /**
-     * A class like the one with same name in tria.cc. See there for more
-     * information.
-     */
-    struct Implementation
-    {
-      /**
-       * In several applications of DoFAccessor::get_dof_values(), we want to
-       * extract some indices without having to allocate memory. We do this by
-       * setting a boost small_vector with 27 elements on the stack, and only
-       * allocate when we exceed 27. The number 27 is heuristic and allows up
-       * to quadratic shape functions on scalar problems in 3d, or linear
-       * shape functions on systems (elasticity).
-       */
-      using dof_index_vector_type =
-        boost::container::small_vector<dealii::types::global_dof_index, 27>;
-
-      /**
-       * Process the @p local_index-th degree of freedom corresponding to the
-       * finite element specified by @p fe_index on the vertex with global
-       * number @p vertex_index to @p global_index.
-       *
-       * The template argument `structdim` indicates the
-       * dimensionality of the object on which we seek to know the DoF
-       * index. For example, if `structdim==0`, then we are looking to
-       * get a DoF index on a vertex of the indicated cell.
-       */
-      template <int dim,
-                int spacedim,
-                int structdim,
-                typename GlobalIndexType,
-                typename DoFPProcessor>
-      static void
-      process_dof_index(const DoFHandler<dim, spacedim> &dof_handler,
-                        const unsigned int               obj_level,
-                        const unsigned int               obj_index,
-                        const types::fe_index            fe_index,
-                        const unsigned int               local_index,
-                        const std::integral_constant<int, structdim> &,
-                        GlobalIndexType     &global_index,
-                        const DoFPProcessor &process)
-      {
-        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
-
-        // 1) no hp used -> fe_index == 0
-        if (dof_handler.hp_capability_enabled == false)
-          {
-            AssertDimension(fe_index,
-                            (DoFHandler<dim, spacedim>::default_fe_index));
-
-            process(
-              dof_handler.object_dof_indices
-                [obj_level][structdim]
-                [dof_handler.object_dof_ptr[obj_level][structdim][obj_index] +
-                 local_index],
-              global_index);
-
-            return;
-          }
-
-        // 2) cell and hp is used -> there is only one fe_index
-        if (structdim == dim)
-          {
-            process(
-              dof_handler.object_dof_indices
-                [obj_level][structdim]
-                [dof_handler.object_dof_ptr[obj_level][structdim][obj_index] +
-                 local_index],
-              global_index);
-            return;
-          }
-
-        // 3) general entity and hp is used
-        AssertIndexRange(obj_level, dof_handler.object_dof_indices.size());
-        AssertIndexRange(structdim,
-                         dof_handler.object_dof_indices[obj_level].size());
-
-        Assert(dof_handler.hp_capability_enabled, ExcInternalError());
-
-        AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
-        AssertIndexRange(obj_index,
-                         dof_handler.hp_object_fe_ptr[structdim].size());
-
-        const auto ptr =
-          std::find(dof_handler.hp_object_fe_indices[structdim].begin() +
-                      dof_handler.hp_object_fe_ptr[structdim][obj_index],
-                    dof_handler.hp_object_fe_indices[structdim].begin() +
-                      dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
-                    fe_index);
-
-        Assert(ptr != dof_handler.hp_object_fe_indices[structdim].begin() +
-                        dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
-               ExcMessage(
-                 "You are requesting an active FE index that is not assigned "
-                 "to any of the cells connected to this entity."));
-
-        const types::fe_index fe_index_ =
-          std::distance(dof_handler.hp_object_fe_indices[structdim].begin() +
-                          dof_handler.hp_object_fe_ptr[structdim][obj_index],
-                        ptr);
-
-        AssertIndexRange(
-          dof_handler.hp_capability_enabled ?
-            (dof_handler.hp_object_fe_ptr[structdim][obj_index] + fe_index_) :
-            obj_index,
-          dof_handler.object_dof_ptr[obj_level][structdim].size());
-
-        AssertIndexRange(
-          dof_handler.object_dof_ptr
-              [obj_level][structdim]
-              [dof_handler.hp_capability_enabled ?
-                 (dof_handler.hp_object_fe_ptr[structdim][obj_index] +
-                  fe_index_) :
-                 obj_index] +
-            local_index,
-          dof_handler.object_dof_indices[obj_level][structdim].size());
-
-        process(dof_handler.object_dof_indices
-                  [obj_level][structdim]
-                  [dof_handler.object_dof_ptr
-                     [obj_level][structdim]
-                     [dof_handler.hp_capability_enabled ?
-                        (dof_handler.hp_object_fe_ptr[structdim][obj_index] +
-                         fe_index_) :
-                        obj_index] +
-                   local_index],
-                global_index);
-      }
-
-      /**
-       * Determine start index and number of dofs of object in global data
-       * structure.
-       *
-       * The template argument `structdim` indicates the
-       * dimensionality of the object on which we seek to know the DoF
-       * index. For example, if `structdim==0`, then we are looking to
-       * get a DoF index on a vertex of the indicated cell.
-       */
-      template <int dim, int spacedim, int structdim>
-      static std::pair<unsigned int, unsigned int>
-      process_object_range(const DoFHandler<dim, spacedim> &dof_handler,
-                           const unsigned int               obj_level,
-                           const unsigned int               obj_index,
-                           const types::fe_index            fe_index,
-                           const std::integral_constant<int, structdim> &)
-      {
-        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
-
-        // determine range of dofs in global data structure
-        // 1) cell
-        if (structdim == dim)
-          {
-            const unsigned int ptr_0 =
-              dof_handler.object_dof_ptr[obj_level][structdim][obj_index];
-            const unsigned int length =
-              dof_handler.get_fe(fe_index).template n_dofs_per_object<dim>(0);
-
-            return {ptr_0, length};
-          }
-
-        // 2) hp is not used -> fe_index == 0
-        if (dof_handler.hp_capability_enabled == false)
-          {
-            AssertDimension(fe_index,
-                            (DoFHandler<dim, spacedim>::default_fe_index));
-
-            const unsigned int ptr_0 =
-              dof_handler.object_dof_ptr[obj_level][structdim][obj_index];
-            const unsigned int length =
-              dof_handler.object_dof_ptr[obj_level][structdim][obj_index + 1] -
-              ptr_0;
-
-            return {ptr_0, length};
-          }
-
-        // 3) hp is used
-        AssertIndexRange(obj_level, dof_handler.object_dof_indices.size());
-        AssertIndexRange(structdim,
-                         dof_handler.object_dof_indices[obj_level].size());
-
-        AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
-        AssertIndexRange(obj_index,
-                         dof_handler.hp_object_fe_ptr[structdim].size());
-
-        const auto fe_index_local_ptr =
-          std::find(dof_handler.hp_object_fe_indices[structdim].begin() +
-                      dof_handler.hp_object_fe_ptr[structdim][obj_index],
-                    dof_handler.hp_object_fe_indices[structdim].begin() +
-                      dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
-                    fe_index);
-
-        Assert(fe_index_local_ptr !=
-                 dof_handler.hp_object_fe_indices[structdim].begin() +
-                   dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
-               ExcMessage(
-                 "You tried to call a function accessing DoF indices, but "
-                 "they appear not be available (yet) or inconsistent. "
-                 "Did you call distribute_dofs() first? Alternatively, if "
-                 "you are using different elements on different cells (i.e., "
-                 "you are using the hp capabilities of deal.II), did you "
-                 "change the active_fe_index of a cell since you last "
-                 "called distribute_dofs()?"));
-
-        const types::fe_index fe_index_local =
-          std::distance(dof_handler.hp_object_fe_indices[structdim].begin() +
-                          dof_handler.hp_object_fe_ptr[structdim][obj_index],
-                        fe_index_local_ptr);
-
-        AssertIndexRange(
-          dof_handler.hp_object_fe_ptr[structdim][obj_index] + fe_index_local,
-          dof_handler.object_dof_ptr[obj_level][structdim].size());
-
-        const unsigned int ptr_0 =
-          dof_handler
-            .object_dof_ptr[obj_level][structdim]
-                           [dof_handler.hp_object_fe_ptr[structdim][obj_index] +
-                            fe_index_local];
-        const unsigned int ptr_1 =
-          dof_handler
-            .object_dof_ptr[obj_level][structdim]
-                           [dof_handler.hp_object_fe_ptr[structdim][obj_index] +
-                            fe_index_local + 1];
-
-        return {ptr_0, ptr_1 - ptr_0};
-      }
-
-      template <int dim, int spacedim, int structdim, bool level_dof_access>
-      static std::pair<unsigned int, unsigned int>
-      process_object_range(
-        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
-                              accessor,
-        const types::fe_index fe_index)
-      {
-        return process_object_range(accessor.get_dof_handler(),
-                                    accessor.level(),
-                                    accessor.index(),
-                                    fe_index,
-                                    std::integral_constant<int, structdim>());
-      }
-
-      template <int dim, int spacedim, int structdim>
-      static std::pair<unsigned int, unsigned int>
-      process_object_range(dealii::DoFInvalidAccessor<structdim, dim, spacedim>,
-                           const unsigned int)
-      {
-        DEAL_II_ASSERT_UNREACHABLE();
-
-        return {0, 0};
-      }
-
-
-
-      /**
-       * Process all dofs of an object.
-       *
-       * The template argument `structdim` indicates the
-       * dimensionality of the object on which we seek to know the DoF
-       * index. For example, if `structdim==0`, then we are looking to
-       * get a DoF index on a vertex of the indicated cell.
-       */
-      template <int dim,
-                int spacedim,
-                int structdim,
-                typename DoFProcessor,
-                typename DoFMapping>
-      static DEAL_II_ALWAYS_INLINE void
-      process_object(const DoFHandler<dim, spacedim>              &dof_handler,
-                     const unsigned int                            obj_level,
-                     const unsigned int                            obj_index,
-                     const types::fe_index                         fe_index,
-                     const DoFMapping                             &mapping,
-                     const std::integral_constant<int, structdim> &dd,
-                     types::global_dof_index *&dof_indices_ptr,
-                     const DoFProcessor       &process)
-      {
-        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
-
-        // determine range of dofs in global data structure
-        const auto range =
-          process_object_range(dof_handler, obj_level, obj_index, fe_index, dd);
-        if (range.second == 0)
-          return;
-
-        std::vector<types::global_dof_index> &object_dof_indices =
-          dof_handler
-            .object_dof_indices[structdim < dim ? 0 : obj_level][structdim];
-        AssertIndexRange(range.first, object_dof_indices.size());
-        types::global_dof_index *DEAL_II_RESTRICT stored_indices =
-          object_dof_indices.data() + range.first;
-
-        // process dofs
-        for (unsigned int i = 0; i < range.second; ++i)
-          {
-            process(
-              stored_indices[(structdim == 0 || structdim == dim) ? i :
-                                                                    mapping(i)],
-              dof_indices_ptr);
-            if (dof_indices_ptr != nullptr)
-              ++dof_indices_ptr;
-          }
-      }
-
-
-
-      /**
-       * Set the @p local_index-th degree of freedom corresponding to the
-       * finite element specified by @p fe_index on the vertex with global
-       * number @p vertex_index to @p global_index.
-       *
-       * The template argument `structdim` indicates the
-       * dimensionality of the object on which we seek to know the DoF
-       * index. For example, if `structdim==0`, then we are looking to
-       * get a DoF index on a vertex of the indicated cell.
-       */
-      template <int dim, int spacedim, int structdim>
-      static void
-      set_dof_index(const DoFHandler<dim, spacedim>              &dof_handler,
-                    const unsigned int                            obj_level,
-                    const unsigned int                            obj_index,
-                    const types::fe_index                         fe_index,
-                    const unsigned int                            local_index,
-                    const std::integral_constant<int, structdim> &dd,
-                    const types::global_dof_index                 global_index)
-      {
-        process_dof_index(dof_handler,
-                          obj_level,
-                          obj_index,
-                          fe_index,
-                          local_index,
-                          dd,
-                          global_index,
-                          [](auto &ptr, const auto &value) { ptr = value; });
-      }
-
-
-      /**
-       * Get the @p local_index-th degree of freedom corresponding to the
-       * finite element specified by @p fe_index on the vertex with global
-       * number @p vertex_index to @p global_index.
-       *
-       * The template argument `structdim` indicates the
-       * dimensionality of the object on which we seek to know the DoF
-       * index. For example, if `structdim==0`, then we are looking to
-       * get a DoF index on a vertex of the indicated cell.
-       */
-      template <int dim, int spacedim, int structdim>
-      static types::global_dof_index
-      get_dof_index(const DoFHandler<dim, spacedim>              &dof_handler,
-                    const unsigned int                            obj_level,
-                    const unsigned int                            obj_index,
-                    const types::fe_index                         fe_index,
-                    const unsigned int                            local_index,
-                    const std::integral_constant<int, structdim> &dd)
-      {
-        types::global_dof_index global_index;
-        process_dof_index(dof_handler,
-                          obj_level,
-                          obj_index,
-                          fe_index,
-                          local_index,
-                          dd,
-                          global_index,
-                          [](const auto &ptr, auto &value) { value = ptr; });
-        return global_index;
-      }
-
-
-      template <int dim, int spacedim>
-      static types::global_dof_index &
-      mg_vertex_dof_index(DoFHandler<dim, spacedim> &dof_handler,
-                          const int                  level,
-                          const unsigned int         vertex_index,
-                          const unsigned int         i)
-      {
-        Assert(dof_handler.hp_capability_enabled == false,
-               ExcMessage(
-                 "DoFHandler in hp-mode does not implement multilevel DoFs."));
-
-        return dof_handler.mg_vertex_dofs[vertex_index].access_index(
-          level, i, dof_handler.get_fe().n_dofs_per_vertex());
-      }
-
-
-
-      /**
-       * Return the number of different finite elements that are active on a
-       * given object such as a vertex, line, or cell.
-       *
-       * The template argument `structdim` indicates the
-       * dimensionality of the object on which we seek to know the DoF
-       * index. For example, if `structdim==0`, then we are looking to
-       * get a DoF index on a vertex of the indicated cell.
-       */
-      template <int dim, int spacedim, int structdim>
-      static unsigned int
-      n_active_fe_indices(const DoFHandler<dim, spacedim> &dof_handler,
-                          const unsigned int               obj_level,
-                          const unsigned int               obj_index,
-                          const std::integral_constant<int, structdim> &)
-      {
-        (void)obj_level;
-
-        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
-
-        // 1) no hp used -> fe_index == 0
-        if (dof_handler.hp_capability_enabled == false)
-          return 1;
-
-        // 2) cell and hp is used -> there is only one fe_index
-        if (structdim == dim)
-          return 1;
-
-        // 3) general entity and hp is used
-        AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
-        AssertIndexRange(obj_index + 1,
-                         dof_handler.hp_object_fe_ptr[structdim].size());
-
-        return dof_handler.hp_object_fe_ptr[structdim][obj_index + 1] -
-               dof_handler.hp_object_fe_ptr[structdim][obj_index];
-      }
-
-
-
-      /**
-       * Return the FE index of the n-th finite element active on a given
-       * object such as a vertex, line, or cell.
-       *
-       * The template argument `structdim` indicates the
-       * dimensionality of the object on which we seek to know the DoF
-       * index. For example, if `structdim==0`, then we are looking to
-       * get a DoF index on a vertex of the indicated cell.
-       */
-      template <int dim, int spacedim, int structdim>
-      static types::fe_index
-      nth_active_fe_index(const DoFHandler<dim, spacedim> &dof_handler,
-                          const unsigned int               obj_level,
-                          const unsigned int               obj_index,
-                          const unsigned int               local_index,
-                          const std::integral_constant<int, structdim> &)
-      {
-        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
-
-        // for cells only one active FE index available
-        Assert(((structdim == dim) &&
-                (local_index != DoFHandler<dim, spacedim>::default_fe_index)) ==
-                 false,
-               ExcNotImplemented());
-
-        // 1) no hp used -> fe_index == 0
-        if (dof_handler.hp_capability_enabled == false)
-          return DoFHandler<dim, spacedim>::default_fe_index;
-
-        // 2) cell and hp is used -> there is only one fe_index
-        if (structdim == dim)
-          return dof_handler.hp_cell_active_fe_indices[obj_level][obj_index];
-
-        // 3) general entity and hp is used
-        AssertIndexRange(structdim, dof_handler.hp_object_fe_indices.size());
-        AssertIndexRange(structdim, dof_handler.hp_object_fe_ptr.size());
-        AssertIndexRange(obj_index,
-                         dof_handler.hp_object_fe_ptr[structdim].size());
-        AssertIndexRange(dof_handler.hp_object_fe_ptr[structdim][obj_index] +
-                           local_index,
-                         dof_handler.hp_object_fe_indices[structdim].size());
-
-        return dof_handler.hp_object_fe_indices
-          [structdim]
-          [dof_handler.hp_object_fe_ptr[structdim][obj_index] + local_index];
-      }
-
-
-
-      /**
-       * Returns all active FE indices on a given object such as a
-       * vertex, line, or cell.
-       *
-       * The size of the returned set equals the number of finite elements that
-       * are active on this vertex.
-       *
-       * The template argument `structdim` indicates the
-       * dimensionality of the object on which we seek to know the DoF
-       * index. For example, if `structdim==0`, then we are looking to
-       * get a DoF index on a vertex of the indicated cell.
-       */
-      template <int dim, int spacedim, int structdim>
-      static std::set<types::fe_index>
-      get_active_fe_indices(const DoFHandler<dim, spacedim> &dof_handler,
-                            const unsigned int               obj_level,
-                            const unsigned int               obj_index,
-                            const std::integral_constant<int, structdim> &t)
-      {
-        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
-
-        // 1) no hp used -> fe_index == 0
-        if (dof_handler.hp_capability_enabled == false)
-          return {DoFHandler<dim, spacedim>::default_fe_index};
-
-        // 2) cell and hp is used -> there is only one fe_index
-        if (structdim == dim)
-          return {dof_handler.hp_cell_active_fe_indices[obj_level][obj_index]};
-
-        // 3) general entity and hp is used
-        std::set<types::fe_index> active_fe_indices;
-        for (unsigned int i = 0;
-             i < n_active_fe_indices(dof_handler, obj_level, obj_index, t);
-             ++i)
-          active_fe_indices.insert(
-            nth_active_fe_index(dof_handler, obj_level, obj_index, i, t));
-        return active_fe_indices;
-      }
-
-
-
-      template <int dim, int spacedim, int structdim>
-      static bool
-      fe_index_is_active(const DoFHandler<dim, spacedim> &dof_handler,
-                         const unsigned int               obj_level,
-                         const unsigned int               obj_index,
-                         const types::fe_index            fe_index,
-                         const std::integral_constant<int, structdim> &)
-      {
-        Assert(structdim == dim || obj_level == 0, ExcNotImplemented());
-
-        // 1) no hp used -> fe_index == 0
-        if (dof_handler.hp_capability_enabled == false)
-          return (fe_index == DoFHandler<dim, spacedim>::default_fe_index);
-
-        // 2) cell and hp is used -> there is only one fe_index
-        if (structdim == dim)
-          return dof_handler.hp_cell_active_fe_indices[obj_level][obj_index] ==
-                 fe_index;
-
-        // 3) general entity and hp is used
-        return std::find(
-                 dof_handler.hp_object_fe_indices[structdim].begin() +
-                   dof_handler.hp_object_fe_ptr[structdim][obj_index],
-                 dof_handler.hp_object_fe_indices[structdim].begin() +
-                   dof_handler.hp_object_fe_ptr[structdim][obj_index + 1],
-                 fe_index) !=
-               (dof_handler.hp_object_fe_indices[structdim].begin() +
-                dof_handler.hp_object_fe_ptr[structdim][obj_index + 1]);
-      }
-
-
-
-      template <typename InputVector, typename ForwardIterator>
-      static void
-      extract_subvector_to(const InputVector             &values,
-                           const types::global_dof_index *cache,
-                           const types::global_dof_index *cache_end,
-                           ForwardIterator                local_values_begin)
-      {
-        values.extract_subvector_to(cache, cache_end, local_values_begin);
-      }
-
-
-
-#ifdef DEAL_II_WITH_TRILINOS
-      static std::vector<unsigned int>
-      sort_indices(const types::global_dof_index *v_begin,
-                   const types::global_dof_index *v_end)
-      {
-        // initialize original index locations
-        std::vector<unsigned int> idx(v_end - v_begin);
-        std::iota(idx.begin(), idx.end(), 0u);
-
-        // sort indices based on comparing values in v
-        std::sort(idx.begin(),
-                  idx.end(),
-                  [&v_begin](unsigned int i1, unsigned int i2) {
-                    return *(v_begin + i1) < *(v_begin + i2);
-                  });
-
-        return idx;
-      }
-
-
-
-#  ifdef DEAL_II_TRILINOS_WITH_TPETRA
-      template <typename ForwardIterator, typename Number, typename MemorySpace>
-      static void
-      extract_subvector_to(
-        const LinearAlgebra::TpetraWrappers::Vector<Number, MemorySpace>
-                                      &values,
-        const types::global_dof_index *cache_begin,
-        const types::global_dof_index *cache_end,
-        ForwardIterator                local_values_begin)
-      {
-        std::vector<unsigned int> sorted_indices_pos =
-          sort_indices(cache_begin, cache_end);
-        const unsigned int cache_size = cache_end - cache_begin;
-        std::vector<types::global_dof_index> cache_indices(cache_size);
-        for (unsigned int i = 0; i < cache_size; ++i)
-          cache_indices[i] = *(cache_begin + sorted_indices_pos[i]);
-
-        IndexSet index_set(cache_indices.back() + 1);
-        index_set.add_indices(cache_indices.begin(), cache_indices.end());
-        index_set.compress();
-        LinearAlgebra::ReadWriteVector<Number> read_write_vector(index_set);
-        read_write_vector.import_elements(values, VectorOperation::insert);
-
-        // Copy the elements from read_write_vector and reorder them.
-        for (unsigned int i = 0; i < cache_size; ++i, ++local_values_begin)
-          *local_values_begin = read_write_vector[sorted_indices_pos[i]];
-      }
-#  endif
-
-
-
-      template <typename ForwardIterator>
-      static void
-      extract_subvector_to(const LinearAlgebra::EpetraWrappers::Vector &values,
-                           const types::global_dof_index *cache_begin,
-                           const types::global_dof_index *cache_end,
-                           ForwardIterator                local_values_begin)
-      {
-        std::vector<unsigned int> sorted_indices_pos =
-          sort_indices(cache_begin, cache_end);
-        const unsigned int cache_size = cache_end - cache_begin;
-        std::vector<types::global_dof_index> cache_indices(cache_size);
-        for (unsigned int i = 0; i < cache_size; ++i)
-          cache_indices[i] = *(cache_begin + sorted_indices_pos[i]);
-
-        IndexSet index_set(cache_indices.back() + 1);
-        index_set.add_indices(cache_indices.begin(), cache_indices.end());
-        index_set.compress();
-        LinearAlgebra::ReadWriteVector<double> read_write_vector(index_set);
-        read_write_vector.import_elements(values, VectorOperation::insert);
-
-        // Copy the elements from read_write_vector and reorder them.
-        for (unsigned int i = 0; i < cache_size; ++i, ++local_values_begin)
-          *local_values_begin = read_write_vector[sorted_indices_pos[i]];
-      }
-#endif
-
-      /**
-       * Loop over all degrees of freedom of the object described by the
-       * provided @p accessor and @p fe_index and count them.
-       */
-      template <int dim, int spacedim, bool level_dof_access, int structdim>
-      static unsigned int
-      n_dof_indices(
-        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
-                             &accessor,
-        const types::fe_index fe_index_,
-        const bool            count_level_dofs)
-      {
-        // note: we cannot rely on the template parameter level_dof_access here,
-        // since the function get_mg_dof_indices()/set_mg_dof_indices() can be
-        // called even if level_dof_access==false.
-        if (count_level_dofs)
-          {
-            const auto &fe = accessor.get_fe(fe_index_);
-
-            const unsigned int                                   //
-              dofs_per_vertex = fe.n_dofs_per_vertex(),          //
-              dofs_per_line   = fe.n_dofs_per_line(),            //
-              dofs_per_quad   = fe.n_dofs_per_quad(0 /*dummy*/), //
-              dofs_per_hex    = fe.n_dofs_per_hex();             //
-
-            unsigned int index = 0;
-
-            // 1) VERTEX dofs
-            index += dofs_per_vertex * accessor.n_vertices();
-
-            // 2) LINE dofs
-            if (structdim == 2 || structdim == 3)
-              index += dofs_per_line * accessor.n_lines();
-
-            // 3) FACE dofs
-            if (structdim == 3)
-              index += dofs_per_quad * accessor.n_faces();
-
-            // 4) INNER dofs
-            const unsigned int interior_dofs =
-              structdim == 1 ? dofs_per_line :
-                               (structdim == 2 ? dofs_per_quad : dofs_per_hex);
-
-            index += interior_dofs;
-
-            return index;
-          }
-        else
-          {
-            const auto fe_index =
-              internal::DoFAccessorImplementation::get_fe_index_or_default(
-                accessor, fe_index_);
-
-            unsigned int index = 0;
-
-            // 1) VERTEX dofs
-            for (const auto vertex : accessor.vertex_indices())
-              index += process_object_range(accessor.get_dof_handler(),
-                                            0,
-                                            accessor.vertex_index(vertex),
-                                            fe_index,
-                                            std::integral_constant<int, 0>())
-                         .second;
-
-            // 2) LINE dofs
-            if (structdim == 2 || structdim == 3)
-              for (const auto line : accessor.line_indices())
-                index +=
-                  process_object_range(*accessor.line(line), fe_index).second;
-
-            // 3) FACE dofs
-            if (structdim == 3)
-              for (const auto face : accessor.face_indices())
-                index +=
-                  process_object_range(*accessor.quad(face), fe_index).second;
-
-            // 4) INNER dofs
-            index += process_object_range(accessor, fe_index).second;
-
-            return index;
-          }
-      }
-
-
-
-      // The next few internal helper functions are needed to support various
-      // DoFIndicesType kinds, e.g. actual vectors of DoFIndices or empty
-      // types that we use when we only want to work on the internally stored
-      // DoFs and never extract any number.
-      template <typename ArrayType>
-      static unsigned int
-      get_array_length(const ArrayType &array)
-      {
-        return array.size();
-      }
-
-      static unsigned int
-      get_array_length(const std::tuple<> &)
-      {
-        return 0;
-      }
-
-      template <typename ArrayType>
-      static types::global_dof_index *
-      get_array_ptr(const ArrayType &array)
-      {
-        return const_cast<types::global_dof_index *>(array.data());
-      }
-
-      static types::global_dof_index *
-      get_array_ptr(const std::tuple<> &)
-      {
-        return nullptr;
-      }
-
-
-
-      /**
-       * Loop over all degrees of freedom of the object described by the
-       * provided @p accessor and @p fe_index and perform the static functions
-       * provided by DoFOperation (set/get) on these.
-       */
-      template <int  dim,
-                int  spacedim,
-                bool level_dof_access,
-                int  structdim,
-                typename DoFIndicesType,
-                typename DoFOperation,
-                typename DoFProcessor>
-      static void
-      process_dof_indices(
-        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
-                             &accessor,
-        const DoFIndicesType &const_dof_indices,
-        const types::fe_index fe_index_,
-        const DoFOperation   &dof_operation,
-        const DoFProcessor   &dof_processor,
-        const bool            count_level_dofs)
-      {
-        const types::fe_index fe_index =
-          internal::DoFAccessorImplementation::get_fe_index_or_default(
-            accessor, fe_index_);
-
-        // we cannot rely on the template parameter level_dof_access here, since
-        // the function get_mg_dof_indices()/set_mg_dof_indices() can be called
-        // even if level_dof_access==false.
-        (void)count_level_dofs;
-
-        const auto &fe = accessor.get_fe(fe_index);
-
-        // we want to pass in rvalue 'std::tuple<>' types as `DoFIndicesType`,
-        // but we need non-const references for std::vector<> types, so get in
-        // a const reference here and immediately cast the constness away -
-        // note that any use of the dereferenced invalid type will result in a
-        // segfault
-        types::global_dof_index *dof_indices_ptr =
-          get_array_ptr(const_dof_indices);
-        types::global_dof_index *end_dof_indices =
-          dof_indices_ptr + get_array_length(const_dof_indices);
-
-        // 1) VERTEX dofs, only step into the functions if we actually have
-        // DoFs on them
-        if (fe.n_dofs_per_vertex() > 0)
-          for (const auto vertex : accessor.vertex_indices())
-            dof_operation.process_vertex_dofs(*accessor.dof_handler,
-                                              accessor.vertex_index(vertex),
-                                              fe_index,
-                                              dof_indices_ptr,
-                                              dof_processor);
-
-        // 2) copy dof numbers from the LINE. for lines with the wrong
-        // orientation (which might occur in 3d), we have already made sure that
-        // we're ok by picking the correct vertices (this happens automatically
-        // in the vertex() function). however, if the line is in wrong
-        // orientation, we look at it in flipped orientation and we will have to
-        // adjust the shape function indices that we see to correspond to the
-        // correct (face/cell-local) ordering.
-        if ((structdim == 2 || structdim == 3) && fe.n_dofs_per_line() > 0)
-          {
-            const auto line_indices = internal::TriaAccessorImplementation::
-              Implementation::get_line_indices_of_cell(accessor);
-            const auto line_orientations =
-              internal::TriaAccessorImplementation::Implementation::
-                get_line_orientations_of_cell(accessor);
-
-            for (const auto line : accessor.line_indices())
-              {
-                const bool line_orientation = line_orientations[line];
-                if (line_orientation)
-                  dof_operation.process_dofs(
-                    accessor.get_dof_handler(),
-                    0,
-                    line_indices[line],
-                    fe_index,
-                    [](const auto d) { return d; },
-                    std::integral_constant<int, 1>(),
-                    dof_indices_ptr,
-                    dof_processor);
-                else
-                  dof_operation.process_dofs(
-                    accessor.get_dof_handler(),
-                    0,
-                    line_indices[line],
-                    fe_index,
-                    [&fe, line_orientation](const auto d) {
-                      return fe.adjust_line_dof_index_for_line_orientation(
-                        d, line_orientation);
-                    },
-                    std::integral_constant<int, 1>(),
-                    dof_indices_ptr,
-                    dof_processor);
-              }
-          }
-
-        // 3) copy dof numbers from the FACE. for faces with the wrong
-        // orientation, we have already made sure that we're ok by picking the
-        // correct lines and vertices (this happens automatically in the line()
-        // and vertex() functions). however, if the face is in wrong
-        // orientation, we look at it in flipped orientation and we will have to
-        // adjust the shape function indices that we see to correspond to the
-        // correct (cell-local) ordering. The same applies, if the face_rotation
-        // or face_orientation is non-standard
-        if (structdim == 3 && fe.max_dofs_per_quad() > 0)
-          for (const auto face_no : accessor.face_indices())
-            {
-              const auto combined_orientation =
-                accessor.combined_face_orientation(face_no);
-              const unsigned int quad_index = accessor.quad_index(face_no);
-              if (combined_orientation ==
-                  ReferenceCell::default_combined_face_orientation())
-                dof_operation.process_dofs(
-                  accessor.get_dof_handler(),
-                  0,
-                  quad_index,
-                  fe_index,
-                  [](const auto d) { return d; },
-                  std::integral_constant<int, 2>(),
-                  dof_indices_ptr,
-                  dof_processor);
-              else
-                dof_operation.process_dofs(
-                  accessor.get_dof_handler(),
-                  0,
-                  quad_index,
-                  fe_index,
-                  [&](const auto d) {
-                    return fe.adjust_quad_dof_index_for_face_orientation(
-                      d, face_no, combined_orientation);
-                  },
-                  std::integral_constant<int, 2>(),
-                  dof_indices_ptr,
-                  dof_processor);
-            }
-
-        // 4) INNER dofs - here we need to make sure that the shortcut to not
-        // run the function does not miss the faces of wedge and pyramid
-        // elements where n_dofs_per_object might not return the largest
-        // possible value
-        if (((dim == 3 && structdim == 2) ?
-               fe.max_dofs_per_quad() :
-               fe.template n_dofs_per_object<structdim>()) > 0)
-          dof_operation.process_dofs(
-            accessor.get_dof_handler(),
-            accessor.level(),
-            accessor.index(),
-            fe_index,
-            [&](const auto d) { return d; },
-            std::integral_constant<int, structdim>(),
-            dof_indices_ptr,
-            dof_processor);
-
-        if (dof_indices_ptr != nullptr)
-          {
-            AssertDimension(n_dof_indices(accessor, fe_index, count_level_dofs),
-                            dof_indices_ptr - get_array_ptr(const_dof_indices));
-          }
-
-        // PM: This is a part that should not be reached since it indicates that
-        // an object (and/or its subobjects) is not active. However,
-        // unfortunately this function is called by
-        // DoFTools::set_periodicity_constraints() indirectly by
-        // get_dof_indices() also for artificial faces to determine if a face
-        // is artificial.
-        types::global_dof_index invalid_index = numbers::invalid_dof_index;
-        for (; dof_indices_ptr < end_dof_indices; ++dof_indices_ptr)
-          dof_processor(invalid_index, dof_indices_ptr);
-      }
-
-
-
-      /**
-       * An internal struct encapsulating the task of getting (vertex)
-       * DoF indices.
-       */
-      template <int dim, int spacedim>
-      struct DoFIndexProcessor
-      {
-        /**
-         * Return vertex DoF indices.
-         */
-        template <typename DoFProcessor>
-        DEAL_II_ALWAYS_INLINE void
-        process_vertex_dofs(DoFHandler<dim, spacedim> &dof_handler,
-                            const unsigned int         vertex_index,
-                            const types::fe_index      fe_index,
-                            types::global_dof_index  *&dof_indices_ptr,
-                            const DoFProcessor        &dof_processor) const
-        {
-          process_object(
-            dof_handler,
-            0,
-            vertex_index,
-            fe_index,
-            [](const auto d) {
-              DEAL_II_ASSERT_UNREACHABLE();
-              return d;
-            },
-            std::integral_constant<int, 0>(),
-            dof_indices_ptr,
-            dof_processor);
-        }
-
-        /**
-         * Return DoF indices for lines, quads, and inner degrees of freedom.
-         */
-        template <int structdim, typename DoFMapping, typename DoFProcessor>
-        DEAL_II_ALWAYS_INLINE void
-        process_dofs(const DoFHandler<dim, spacedim> &dof_handler,
-                     const unsigned int               obj_level,
-                     const unsigned int               obj_index,
-                     const types::fe_index            fe_index,
-                     const DoFMapping                &mapping,
-                     const std::integral_constant<int, structdim>,
-                     types::global_dof_index *&dof_indices_ptr,
-                     const DoFProcessor       &dof_processor) const
-        {
-          process_object(
-            dof_handler,
-            obj_level,
-            obj_index,
-            fe_index,
-            mapping,
-            std::integral_constant<int, std::min(structdim, dim)>(),
-            dof_indices_ptr,
-            dof_processor);
-        }
-      };
-
-
-
-      /**
-       * An internal struct encapsulating the task of getting level (vertex)
-       * DoF indices.
-       */
-      template <int dim, int spacedim>
-      struct MGDoFIndexProcessor
-      {
-        /**
-         * Constructor.
-         */
-        MGDoFIndexProcessor(const unsigned int level)
-          : level(level)
-        {}
-
-        /**
-         * Return vertex DoF indices.
-         */
-        template <typename DoFProcessor>
-        DEAL_II_ALWAYS_INLINE void
-        process_vertex_dofs(DoFHandler<dim, spacedim> &dof_handler,
-                            const unsigned int         vertex_index,
-                            const types::fe_index,
-                            types::global_dof_index *&dof_indices_ptr,
-                            const DoFProcessor       &dof_processor) const
-        {
-          const unsigned int n_indices =
-            dof_handler.get_fe(0).template n_dofs_per_object<0>();
-          types::global_dof_index *stored_indices =
-            &dof_handler.mg_vertex_dofs[vertex_index].access_index(level,
-                                                                   0,
-                                                                   n_indices);
-          for (unsigned int d = 0; d < n_indices; ++d, ++dof_indices_ptr)
-            dof_processor(stored_indices[d], dof_indices_ptr);
-        }
-
-        /**
-         * Return DoF indices for lines, quads, and inner degrees of freedom.
-         */
-        template <int structdim, typename DoFMapping, typename DoFProcessor>
-        DEAL_II_ALWAYS_INLINE void
-        process_dofs(const DoFHandler<dim, spacedim> &dof_handler,
-                     const unsigned int,
-                     const unsigned int    obj_index,
-                     const types::fe_index fe_index,
-                     const DoFMapping     &mapping,
-                     const std::integral_constant<int, structdim>,
-                     types::global_dof_index *&dof_indices_ptr,
-                     const DoFProcessor       &dof_processor) const
-        {
-          const unsigned int n_indices =
-            dof_handler.get_fe(0).template n_dofs_per_object<structdim>();
-          types::global_dof_index *stored_indices = &get_mg_dof_index(
-            dof_handler,
-            dof_handler.mg_levels[level],
-            dof_handler.mg_faces,
-            obj_index,
-            fe_index,
-            0,
-            std::integral_constant<int, std::min(structdim, dim)>());
-          for (unsigned int d = 0; d < n_indices; ++d, ++dof_indices_ptr)
-            dof_processor(stored_indices[structdim < dim ? mapping(d) : d],
-                          dof_indices_ptr);
-        }
-
-      private:
-        const unsigned int level;
-      };
-
-
-
-      template <int dim, int spacedim, bool level_dof_access, int structdim>
-      static void
-      get_dof_indices(
-        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
-                                             &accessor,
-        std::vector<types::global_dof_index> &dof_indices,
-        const types::fe_index                 fe_index)
-      {
-        process_dof_indices(
-          accessor,
-          dof_indices,
-          fe_index,
-          DoFIndexProcessor<dim, spacedim>(),
-          [](auto stored_index, auto dof_ptr) { *dof_ptr = stored_index; },
-          false);
-      }
-
-
-
-      template <int dim, int spacedim, bool level_dof_access, int structdim>
-      static void
-      set_dof_indices(
-        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
-                                                   &accessor,
-        const std::vector<types::global_dof_index> &dof_indices,
-        const types::fe_index                       fe_index)
-      {
-        // Note: this function is as general as `get_dof_indices()`. This
-        // assert is placed here since it is currently only used by the
-        // function DoFCellAccessor::set_dof_indices(), which is called by
-        // internal::DoFHandlerImplementation::Policy::Implementation::distribute_dofs().
-        // In the case of new use cases, this assert can be removed.
-        Assert(
-          dim == structdim,
-          ExcMessage(
-            "This function is intended to be used for DoFCellAccessor, i.e., "
-            "dimension == structdim."));
-
-        process_dof_indices(
-          accessor,
-          dof_indices,
-          fe_index,
-          DoFIndexProcessor<dim, spacedim>(),
-          [](auto &stored_index, auto dof_ptr) { stored_index = *dof_ptr; },
-          false);
-      }
-
-
-
-      template <int dim, int spacedim, bool level_dof_access, int structdim>
-      static void
-      get_mg_dof_indices(
-        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
-                                             &accessor,
-        const int                             level,
-        std::vector<types::global_dof_index> &dof_indices,
-        const types::fe_index                 fe_index)
-      {
-        Assert((fe_index == DoFHandler<dim, spacedim>::default_fe_index),
-               ExcMessage("MG DoF indices cannot be queried in hp case"));
-        process_dof_indices(
-          accessor,
-          dof_indices,
-          fe_index,
-          MGDoFIndexProcessor<dim, spacedim>(level),
-          [](auto stored_index, auto dof_ptr) { *dof_ptr = stored_index; },
-          true);
-      }
-
-
-
-      template <int dim, int spacedim, bool level_dof_access, int structdim>
-      static void
-      set_mg_dof_indices(
-        const dealii::DoFAccessor<structdim, dim, spacedim, level_dof_access>
-                                                   &accessor,
-        const int                                   level,
-        const std::vector<types::global_dof_index> &dof_indices,
-        const types::fe_index                       fe_index)
-      {
-        Assert((fe_index == DoFHandler<dim, spacedim>::default_fe_index),
-               ExcMessage("MG DoF indices cannot be queried in hp case"));
-
-        // Note: this function is as general as `get_mg_dof_indices()`. This
-        // assert is placed here since it is currently only used by the
-        // function DoFCellAccessor::set_mg_dof_indices(), which is called by
-        // internal::DoFHandlerImplementation::Policy::Implementation::distribute_mg_dofs().
-        // In the case of new use cases, this assert can be removed.
-        Assert(dim == structdim,
-               ExcMessage("This function is intended to be used for "
-                          "DoFCellAccessor, i.e., dimension == structdim."));
-
-        process_dof_indices(
-          accessor,
-          dof_indices,
-          fe_index,
-          MGDoFIndexProcessor<dim, spacedim>(level),
-          [](auto &stored_index, auto dof_ptr) { stored_index = *dof_ptr; },
-          true);
-      }
-
-
-
-      template <int dim, int spacedim>
-      static types::global_dof_index &
-      get_mg_dof_index(
-        const DoFHandler<dim, spacedim> &dof_handler,
-        const std::unique_ptr<internal::DoFHandlerImplementation::DoFLevel<dim>>
-          &mg_level,
-        const std::unique_ptr<internal::DoFHandlerImplementation::DoFFaces<dim>>
-          &,
-        const unsigned int    obj_index,
-        const types::fe_index fe_index,
-        const unsigned int    local_index,
-        const std::integral_constant<int, dim>)
-      {
-        Assert(dof_handler.hp_capability_enabled == false,
-               (typename DoFHandler<dim, spacedim>::ExcNotImplementedWithHP()));
-
-        return mg_level->dof_object.access_dof_index(
-          static_cast<const DoFHandler<dim, spacedim> &>(dof_handler),
-          obj_index,
-          fe_index,
-          local_index);
-      }
-
-
-
-      template <int dim, int spacedim, std::enable_if_t<(dim > 1), int> = 0>
-      static types::global_dof_index &
-      get_mg_dof_index(
-        const DoFHandler<dim, spacedim> &dof_handler,
-        const std::unique_ptr<internal::DoFHandlerImplementation::DoFLevel<dim>>
-          &,
-        const std::unique_ptr<internal::DoFHandlerImplementation::DoFFaces<dim>>
-                             &mg_faces,
-        const unsigned int    obj_index,
-        const types::fe_index fe_index,
-        const unsigned int    local_index,
-        const std::integral_constant<int, 1>)
-      {
-        return mg_faces->lines.access_dof_index(
-          static_cast<const DoFHandler<dim, spacedim> &>(dof_handler),
-          obj_index,
-          fe_index,
-          local_index);
-      }
-
-
-
-      template <int spacedim>
-      static types::global_dof_index &
-      get_mg_dof_index(
-        const DoFHandler<3, spacedim> &dof_handler,
-        const std::unique_ptr<internal::DoFHandlerImplementation::DoFLevel<3>>
-          &,
-        const std::unique_ptr<internal::DoFHandlerImplementation::DoFFaces<3>>
-                             &mg_faces,
-        const unsigned int    obj_index,
-        const types::fe_index fe_index,
-        const unsigned int    local_index,
-        const std::integral_constant<int, 2>)
-      {
-        Assert(dof_handler.hp_capability_enabled == false,
-               (typename DoFHandler<3, spacedim>::ExcNotImplementedWithHP()));
-        return mg_faces->quads.access_dof_index(
-          static_cast<const DoFHandler<3, spacedim> &>(dof_handler),
-          obj_index,
-          fe_index,
-          local_index);
-      }
-    };
-
-
-
-    template <int dim, int spacedim, bool level_dof_access>
-    void
-    get_cell_dof_indices(
-      const dealii::DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
-      Implementation::dof_index_vector_type &dof_indices,
-      const unsigned int                     fe_index);
-  } // namespace DoFAccessorImplementation
-} // namespace internal
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline types::global_dof_index
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::dof_index(
-  const unsigned int    i,
-  const types::fe_index fe_index_) const
-{
-  const auto fe_index =
-    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
-                                                                 fe_index_);
-
-  // access the respective DoF
-  return dealii::internal::DoFAccessorImplementation::Implementation::
-    get_dof_index(*this->dof_handler,
-                  this->level(),
-                  this->index(),
-                  fe_index,
-                  i,
-                  std::integral_constant<int, structdim>());
-}
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline types::global_dof_index
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::mg_dof_index(
-  const int          level,
-  const unsigned int i) const
-{
-  return internal::DoFAccessorImplementation::Implementation::get_mg_dof_index(
-    *this->dof_handler,
-    this->dof_handler->mg_levels[level],
-    this->dof_handler->mg_faces,
-    this->index(),
-    0,
-    i,
-    std::integral_constant<int, structdim>());
-}
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline void
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::set_dof_index(
-  const unsigned int            i,
-  const types::global_dof_index index,
-  const types::fe_index         fe_index_) const
-{
-  const auto fe_index =
-    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
-                                                                 fe_index_);
-
-  // access the respective DoF
-  dealii::internal::DoFAccessorImplementation::Implementation::set_dof_index(
-    *this->dof_handler,
-    this->level(),
-    this->index(),
-    fe_index,
-    i,
-    std::integral_constant<int, structdim>(),
-    index);
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline void
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::set_mg_dof_index(
-  const int                     level,
-  const unsigned int            i,
-  const types::global_dof_index index) const
-{
-  internal::DoFAccessorImplementation::Implementation::get_mg_dof_index(
-    *this->dof_handler,
-    this->dof_handler->mg_levels[level],
-    this->dof_handler->mg_faces,
-    this->index(),
-    0,
-    i,
-    std::integral_constant<int, structdim>()) = index;
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline unsigned int
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::n_active_fe_indices()
-  const
-{
-  // access the respective DoF
-  return dealii::internal::DoFAccessorImplementation::Implementation::
-    n_active_fe_indices(*this->dof_handler,
-                        this->level(),
-                        this->index(),
-                        std::integral_constant<int, structdim>());
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline types::fe_index
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::nth_active_fe_index(
-  const unsigned int n) const
-{
-  // access the respective DoF
-  return dealii::internal::DoFAccessorImplementation::Implementation::
-    nth_active_fe_index(*this->dof_handler,
-                        this->level(),
-                        this->index(),
-                        n,
-                        std::integral_constant<int, structdim>());
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline std::set<types::fe_index>
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::get_active_fe_indices()
-  const
-{
-  std::set<types::fe_index> active_fe_indices;
-  for (unsigned int i = 0; i < n_active_fe_indices(); ++i)
-    active_fe_indices.insert(nth_active_fe_index(i));
-  return active_fe_indices;
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline bool
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::fe_index_is_active(
-  const types::fe_index fe_index) const
-{
-  // access the respective DoF
-  return dealii::internal::DoFAccessorImplementation::Implementation::
-    fe_index_is_active(*this->dof_handler,
-                       this->level(),
-                       this->index(),
-                       fe_index,
-                       std::integral_constant<int, structdim>());
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline types::global_dof_index
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::vertex_dof_index(
-  const unsigned int    vertex,
-  const unsigned int    i,
-  const types::fe_index fe_index_) const
-{
-  const types::fe_index fe_index =
-    (((this->dof_handler->hp_capability_enabled == false) &&
-      (fe_index_ == numbers::invalid_fe_index)) ?
-       // No hp enabled, and the argument is at its default value -> we
-       // can translate to the default active fe index
-       DoFHandler<dim, spacedim>::default_fe_index :
-       // Otherwise: If anything other than the default is provided by
-       // the caller, then we should take just that. As an exception, if
-       // we are on a cell (rather than a face/edge/vertex), then we know
-       // that there is only one active fe index on this cell and we can
-       // use that:
-       ((dim == structdim) && (fe_index_ == numbers::invalid_fe_index) ?
-          this->nth_active_fe_index(0) :
-          fe_index_));
-
-  return dealii::internal::DoFAccessorImplementation::Implementation::
-    get_dof_index(*this->dof_handler,
-                  0,
-                  this->vertex_index(vertex),
-                  fe_index,
-                  i,
-                  std::integral_constant<int, 0>());
-}
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline types::global_dof_index
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::mg_vertex_dof_index(
-  const int             level,
-  const unsigned int    vertex,
-  const unsigned int    i,
-  const types::fe_index fe_index_) const
-{
-  const auto fe_index =
-    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
-                                                                 fe_index_);
-  (void)fe_index;
-  Assert(this->dof_handler != nullptr, ExcInvalidObject());
-  Assert(this->dof_handler->mg_vertex_dofs.size() > 0,
-         ExcMessage("Multigrid DoF indices can only be accessed after "
-                    "DoFHandler::distribute_mg_dofs() has been called!"));
-  AssertIndexRange(vertex, this->n_vertices());
-  AssertIndexRange(i, this->dof_handler->get_fe(fe_index).n_dofs_per_vertex());
-
-  Assert(dof_handler->hp_capability_enabled == false,
-         ExcMessage(
-           "DoFHandler in hp-mode does not implement multilevel DoFs."));
-
-  return this->dof_handler->mg_vertex_dofs[this->vertex_index(vertex)]
-    .access_index(level, i, this->dof_handler->get_fe().n_dofs_per_vertex());
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline void
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::
-  set_mg_vertex_dof_index(const int                     level,
-                          const unsigned int            vertex,
-                          const unsigned int            i,
-                          const types::global_dof_index index,
-                          const types::fe_index         fe_index_) const
-{
-  const auto fe_index =
-    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
-                                                                 fe_index_);
-  (void)fe_index;
-  Assert(this->dof_handler != nullptr, ExcInvalidObject());
-  AssertIndexRange(vertex, this->n_vertices());
-  AssertIndexRange(i, this->dof_handler->get_fe(fe_index).n_dofs_per_vertex());
-
-  Assert(dof_handler->hp_capability_enabled == false,
-         ExcMessage(
-           "DoFHandler in hp-mode does not implement multilevel DoFs."));
-
-  this->dof_handler->mg_vertex_dofs[this->vertex_index(vertex)].access_index(
-    level, i, this->dof_handler->get_fe().n_dofs_per_vertex()) = index;
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline const FiniteElement<dim, spacedim> &
-DoFAccessor<structdim, dim, spacedim, level_dof_access>::get_fe(
-  const types::fe_index fe_index) const
-{
-  Assert(fe_index_is_active(fe_index) == true,
-         ExcMessage("This function can only be called for active FE indices"));
-
-  return this->dof_handler->get_fe(fe_index);
-}
-
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline typename dealii::internal::DoFHandlerImplementation::
-  Iterators<dim, spacedim, level_dof_access>::line_iterator
-  DoFAccessor<structdim, dim, spacedim, level_dof_access>::line(
-    const unsigned int i) const
-{
-  // if we are asking for a particular line and this object refers to
-  // a line, then the only valid index is i==0 and we should return
-  // *this
-  if (structdim == 1)
-    {
-      Assert(i == 0,
-             ExcMessage("You can only ask for line zero if the "
-                        "current object is a line itself."));
-      return typename dealii::internal::DoFHandlerImplementation::
-        Iterators<dim, spacedim, level_dof_access>::cell_iterator(
-          &this->get_triangulation(),
-          this->level(),
-          this->index(),
-          &this->get_dof_handler());
-    }
-
-  // otherwise we need to be in structdim>=2
-  Assert(structdim > 1, ExcImpossibleInDim(structdim));
-  Assert(dim > 1, ExcImpossibleInDim(dim));
-
-  // checking of 'i' happens in line_index(i)
-  return typename dealii::internal::DoFHandlerImplementation::
-    Iterators<dim, spacedim, level_dof_access>::line_iterator(
-      this->tria,
-      0, // only sub-objects are allowed, which have no level
-      this->line_index(i),
-      this->dof_handler);
-}
-
-
-template <int structdim, int dim, int spacedim, bool level_dof_access>
-inline typename dealii::internal::DoFHandlerImplementation::
-  Iterators<dim, spacedim, level_dof_access>::quad_iterator
-  DoFAccessor<structdim, dim, spacedim, level_dof_access>::quad(
-    const unsigned int i) const
-{
-  // if we are asking for a
-  // particular quad and this object
-  // refers to a quad, then the only
-  // valid index is i==0 and we
-  // should return *this
-  if (structdim == 2)
-    {
-      Assert(i == 0,
-             ExcMessage("You can only ask for quad zero if the "
-                        "current object is a quad itself."));
-      return typename dealii::internal::DoFHandlerImplementation::
-        Iterators<dim, spacedim>::cell_iterator(&this->get_triangulation(),
-                                                this->level(),
-                                                this->index(),
-                                                &this->get_dof_handler());
-    }
-
-  // otherwise we need to be in structdim>=3
-  Assert(structdim > 2, ExcImpossibleInDim(structdim));
-  Assert(dim > 2, ExcImpossibleInDim(dim));
-
-  // checking of 'i' happens in quad_index(i)
-  return typename dealii::internal::DoFHandlerImplementation::
-    Iterators<dim, spacedim, level_dof_access>::quad_iterator(
-      this->tria,
-      0, // only sub-objects are allowed, which have no level
-      this->quad_index(i),
-      this->dof_handler);
-}
-
-
-/*----------------- Functions: DoFAccessor<0,1,spacedim> --------------------*/
-
-
-template <int spacedim, bool level_dof_access>
-inline DoFAccessor<0, 1, spacedim, level_dof_access>::DoFAccessor()
-{
-  Assert(false, ExcInvalidObject());
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline DoFAccessor<0, 1, spacedim, level_dof_access>::DoFAccessor(
-  const Triangulation<1, spacedim>                       *tria,
-  const typename TriaAccessor<0, 1, spacedim>::VertexKind vertex_kind,
-  const unsigned int                                      vertex_index,
-  const DoFHandler<1, spacedim>                          *dof_handler)
-  : BaseClass(tria, vertex_kind, vertex_index)
-  , dof_handler(const_cast<DoFHandler<1, spacedim> *>(dof_handler))
-{}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline DoFAccessor<0, 1, spacedim, level_dof_access>::DoFAccessor(
-  const Triangulation<1, spacedim> *tria,
-  const int                         level,
-  const int                         index,
-  const DoFHandler<1, spacedim>    *dof_handler)
-  // This is the constructor signature for "ordinary" (non-vertex)
-  // accessors and we shouldn't be calling it altogether. But it is also
-  // the constructor that the default-constructor of TriaRawIterator
-  // calls when default-constructing an iterator object. If so, this
-  // happens with level==-2 and index==-2, and this is the only case we
-  // would like to support. We do this by just forwarding to the
-  // other constructor of this class, and then asserting the condition
-  // on level and index.
-  : DoFAccessor<0, 1, spacedim, level_dof_access>(
-      tria,
-      TriaAccessor<0, 1, spacedim>::interior_vertex,
-      0U,
-      dof_handler)
-{
-  (void)level;
-  (void)index;
-  Assert((tria == nullptr) && (level == -2) && (index == -2) &&
-           (dof_handler == nullptr),
-         ExcMessage(
-           "This constructor can not be called for face iterators in 1d, "
-           "except to default-construct iterator objects."));
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-template <int structdim2, int dim2, int spacedim2>
-DoFAccessor<0, 1, spacedim, level_dof_access>::DoFAccessor(
-  const InvalidAccessor<structdim2, dim2, spacedim2> &)
-{
-  Assert(false, ExcInvalidObject());
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
-inline DoFAccessor<0, 1, spacedim, level_dof_access>::DoFAccessor(
-  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &)
-{
-  Assert(false, ExcInvalidObject());
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline void
-DoFAccessor<0, 1, spacedim, level_dof_access>::set_dof_handler(
-  DoFHandler<1, spacedim> *dh)
-{
-  Assert(dh != nullptr, ExcInvalidObject());
-  this->dof_handler = dh;
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline void
-DoFAccessor<0, 1, spacedim, level_dof_access>::set_dof_index(
-  const unsigned int /*i*/,
-  const types::global_dof_index /*index*/,
-  const types::fe_index /*fe_index*/) const
-{
-  DEAL_II_NOT_IMPLEMENTED();
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline const DoFHandler<1, spacedim> &
-DoFAccessor<0, 1, spacedim, level_dof_access>::get_dof_handler() const
-{
-  return *this->dof_handler;
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline void
-DoFAccessor<0, 1, spacedim, level_dof_access>::get_dof_indices(
-  std::vector<types::global_dof_index> &dof_indices,
-  const types::fe_index                 fe_index_) const
-{
-  const auto fe_index =
-    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
-                                                                 fe_index_);
-
-  for (unsigned int i = 0; i < dof_indices.size(); ++i)
-    dof_indices[i] = dealii::internal::DoFAccessorImplementation::
-      Implementation::get_dof_index(*dof_handler,
-                                    0,
-                                    this->global_vertex_index,
-                                    fe_index,
-                                    i,
-                                    std::integral_constant<int, 0>());
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline void
-DoFAccessor<0, 1, spacedim, level_dof_access>::get_mg_dof_indices(
-  const int                             level,
-  std::vector<types::global_dof_index> &dof_indices,
-  const types::fe_index                 fe_index_) const
-{
-  const auto fe_index =
-    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
-                                                                 fe_index_);
-  (void)fe_index;
-  AssertDimension(fe_index, (DoFHandler<1, spacedim>::default_fe_index));
-
-  for (unsigned int i = 0; i < dof_indices.size(); ++i)
-    dof_indices[i] =
-      dealii::internal::DoFAccessorImplementation::Implementation::
-        mg_vertex_dof_index(*dof_handler, level, this->global_vertex_index, i);
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline types::global_dof_index
-DoFAccessor<0, 1, spacedim, level_dof_access>::vertex_dof_index(
-  const unsigned int    vertex,
-  const unsigned int    i,
-  const types::fe_index fe_index_) const
-{
-  const auto fe_index =
-    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
-                                                                 fe_index_);
-
-  (void)vertex;
-  AssertIndexRange(vertex, 1);
-  return dealii::internal::DoFAccessorImplementation::Implementation::
-    get_dof_index(*dof_handler,
-                  0,
-                  this->global_vertex_index,
-                  fe_index,
-                  i,
-                  std::integral_constant<int, 0>());
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline types::global_dof_index
-DoFAccessor<0, 1, spacedim, level_dof_access>::dof_index(
-  const unsigned int    i,
-  const types::fe_index fe_index_) const
-{
-  const auto fe_index =
-    internal::DoFAccessorImplementation::get_fe_index_or_default(*this,
-                                                                 fe_index_);
-
-  return dealii::internal::DoFAccessorImplementation::Implementation::
-    get_dof_index(*this->dof_handler,
-                  0,
-                  this->vertex_index(0),
-                  fe_index,
-                  i,
-                  std::integral_constant<int, 0>());
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline unsigned int
-DoFAccessor<0, 1, spacedim, level_dof_access>::n_active_fe_indices() const
-{
-  return 1;
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline types::fe_index
-DoFAccessor<0, 1, spacedim, level_dof_access>::nth_active_fe_index(
-  const unsigned int /*n*/) const
-{
-  return 0;
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline bool
-DoFAccessor<0, 1, spacedim, level_dof_access>::fe_index_is_active(
-  const types::fe_index /*fe_index*/) const
-{
-  DEAL_II_NOT_IMPLEMENTED();
-  return false;
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline const FiniteElement<1, spacedim> &
-DoFAccessor<0, 1, spacedim, level_dof_access>::get_fe(
-  const types::fe_index fe_index) const
-{
-  Assert(this->dof_handler != nullptr, ExcInvalidObject());
-  return dof_handler->get_fe(fe_index);
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline void
-DoFAccessor<0, 1, spacedim, level_dof_access>::copy_from(
-  const TriaAccessorBase<0, 1, spacedim> &da)
-{
-  Assert(this->dof_handler != nullptr, ExcInvalidObject());
-  BaseClass::copy_from(da);
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-template <bool level_dof_access2>
-inline void
-DoFAccessor<0, 1, spacedim, level_dof_access>::copy_from(
-  const DoFAccessor<0, 1, spacedim, level_dof_access2> &a)
-{
-  BaseClass::copy_from(a);
-  set_dof_handler(a.dof_handler);
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline TriaIterator<DoFAccessor<0, 1, spacedim, level_dof_access>>
-DoFAccessor<0, 1, spacedim, level_dof_access>::child(
-  const unsigned int /*i*/) const
-{
-  return TriaIterator<DoFAccessor<0, 1, spacedim, level_dof_access>>();
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline typename dealii::internal::DoFHandlerImplementation::
-  Iterators<1, spacedim, level_dof_access>::line_iterator
-  DoFAccessor<0, 1, spacedim, level_dof_access>::line(
-    const unsigned int /*c*/) const
-{
-  DEAL_II_NOT_IMPLEMENTED();
-  return typename dealii::internal::DoFHandlerImplementation::
-    Iterators<1, spacedim, level_dof_access>::line_iterator();
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-inline typename dealii::internal::DoFHandlerImplementation::
-  Iterators<1, spacedim, level_dof_access>::quad_iterator
-  DoFAccessor<0, 1, spacedim, level_dof_access>::quad(
-    const unsigned int /*c*/) const
-{
-  DEAL_II_NOT_IMPLEMENTED();
-  return typename dealii::internal::DoFHandlerImplementation::
-    Iterators<1, spacedim, level_dof_access>::quad_iterator();
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
-inline bool
-DoFAccessor<0, 1, spacedim, level_dof_access>::operator==(
-  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &a) const
-{
-  Assert(structdim2 == 0, ExcCantCompareIterators());
-  Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
-  return (BaseClass::operator==(a));
-}
-
-
-
-template <int spacedim, bool level_dof_access>
-template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
-inline bool
-DoFAccessor<0, 1, spacedim, level_dof_access>::operator!=(
-  const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &a) const
-{
-  Assert(structdim2 == 0, ExcCantCompareIterators());
-  Assert(this->dof_handler == a.dof_handler, ExcCantCompareIterators());
-  return (BaseClass::operator!=(a));
-}
-
-
-
-/*------------------------- Functions: DoFCellAccessor -----------------------*/
-
-
-namespace internal
-{
-  namespace DoFCellAccessorImplementation
-  {
-    /**
-     * A class with the same purpose as the similarly named class of the
-     * Triangulation class. See there for more information.
-     */
-    struct Implementation
-    {
-      /**
-       * Do what the active_fe_index function in the parent class is supposed to
-       * do.
-       */
-      template <int dim, int spacedim, bool level_dof_access>
-      static types::fe_index
-      active_fe_index(
-        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor)
-      {
-        if (accessor.dof_handler->hp_capability_enabled == false)
-          return DoFHandler<dim, spacedim>::default_fe_index;
-
-        Assert(accessor.dof_handler != nullptr,
-               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
-        Assert(static_cast<unsigned int>(accessor.level()) <
-                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
-               ExcMessage("DoFHandler not initialized"));
-
-        return accessor.dof_handler
-          ->hp_cell_active_fe_indices[accessor.level()][accessor.present_index];
-      }
-
-
-
-      /**
-       * Do what the set_active_fe_index function in the parent class is
-       * supposed to do.
-       */
-      template <int dim, int spacedim, bool level_dof_access>
-      static void
-      set_active_fe_index(
-        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
-        const types::fe_index                                   i)
-      {
-        if (accessor.dof_handler->hp_capability_enabled == false)
-          {
-            AssertDimension(i, (DoFHandler<dim, spacedim>::default_fe_index));
-            return;
-          }
-
-        Assert(accessor.dof_handler != nullptr,
-               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
-        Assert(static_cast<unsigned int>(accessor.level()) <
-                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
-               ExcMessage("DoFHandler not initialized"));
-        Assert(i != numbers::invalid_fe_index,
-               ExcMessage("Invalid finite element index."));
-
-        accessor.dof_handler
-          ->hp_cell_active_fe_indices[accessor.level()]
-                                     [accessor.present_index] = i;
-      }
-
-
-
-      /**
-       * Do what the future_fe_index function in the parent class is supposed to
-       * do.
-       */
-      template <int dim, int spacedim, bool level_dof_access>
-      static types::fe_index
-      future_fe_index(
-        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor)
-      {
-        if (accessor.dof_handler->hp_capability_enabled == false)
-          return DoFHandler<dim, spacedim>::default_fe_index;
-
-        Assert(accessor.dof_handler != nullptr,
-               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
-        Assert(static_cast<unsigned int>(accessor.level()) <
-                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
-               ExcMessage("DoFHandler not initialized"));
-
-        if (future_fe_index_set(accessor))
-          return accessor.dof_handler
-            ->hp_cell_future_fe_indices[accessor.level()]
-                                       [accessor.present_index];
-        else
-          return accessor.dof_handler
-            ->hp_cell_active_fe_indices[accessor.level()]
-                                       [accessor.present_index];
-      }
-
-
-      /**
-       * Do what the set_future_fe_index function in the parent class is
-       * supposed to do.
-       */
-      template <int dim, int spacedim, bool level_dof_access>
-      static void
-      set_future_fe_index(
-        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor,
-        const types::fe_index                                   i)
-      {
-        if (accessor.dof_handler->hp_capability_enabled == false)
-          {
-            AssertDimension(i, (DoFHandler<dim, spacedim>::default_fe_index));
-            return;
-          }
-
-        Assert(accessor.dof_handler != nullptr,
-               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
-        Assert(static_cast<unsigned int>(accessor.level()) <
-                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
-               ExcMessage("DoFHandler not initialized"));
-        Assert(i != numbers::invalid_fe_index,
-               ExcMessage("Invalid finite element index."));
-
-        accessor.dof_handler
-          ->hp_cell_future_fe_indices[accessor.level()]
-                                     [accessor.present_index] = i;
-      }
-
-
-
-      /**
-       * Do what the future_fe_index_set function in the parent class is
-       * supposed to do.
-       */
-      template <int dim, int spacedim, bool level_dof_access>
-      static bool
-      future_fe_index_set(
-        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor)
-      {
-        if (accessor.dof_handler->hp_capability_enabled == false)
-          return false;
-
-        Assert(accessor.dof_handler != nullptr,
-               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
-        Assert(static_cast<unsigned int>(accessor.level()) <
-                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
-               ExcMessage("DoFHandler not initialized"));
-
-        return accessor.dof_handler
-                 ->hp_cell_future_fe_indices[accessor.level()]
-                                            [accessor.present_index] !=
-               numbers::invalid_fe_index;
-      }
-
-
-
-      /**
-       * Do what the clear_fe_index function in the parent class is supposed to
-       * do.
-       */
-      template <int dim, int spacedim, bool level_dof_access>
-      static void
-      clear_future_fe_index(
-        const DoFCellAccessor<dim, spacedim, level_dof_access> &accessor)
-      {
-        if (accessor.dof_handler->hp_capability_enabled == false)
-          return;
-
-        Assert(accessor.dof_handler != nullptr,
-               (typename std::decay_t<decltype(accessor)>::ExcInvalidObject()));
-        Assert(static_cast<unsigned int>(accessor.level()) <
-                 accessor.dof_handler->hp_cell_future_fe_indices.size(),
-               ExcMessage("DoFHandler not initialized"));
-
-        accessor.dof_handler
-          ->hp_cell_future_fe_indices[accessor.level()]
-                                     [accessor.present_index] =
-          numbers::invalid_fe_index;
-      }
-    };
-  } // namespace DoFCellAccessorImplementation
-} // namespace internal
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  DoFCellAccessor(const Triangulation<dimension_, space_dimension_> *tria,
-                  const int                                          level,
-                  const int                                          index,
-                  const AccessorData                                *local_data)
-  : DoFAccessor<dimension_, dimension_, space_dimension_, level_dof_access>(
-      tria,
-      level,
-      index,
-      local_data)
-{}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <int structdim2, int dim2, int spacedim2>
-inline DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  DoFCellAccessor(const InvalidAccessor<structdim2, dim2, spacedim2> &)
-{
-  Assert(false, typename BaseClass::ExcInvalidObject());
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <int structdim2, int dim2, int spacedim2, bool level_dof_access2>
-inline DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  DoFCellAccessor(
-    const DoFAccessor<structdim2, dim2, spacedim2, level_dof_access2> &other)
-  : BaseClass(other)
-{}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline TriaIterator<
-  DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::neighbor(
-  const unsigned int i) const
-{
-  TriaIterator<DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
-    q(this->tria,
-      this->neighbor_level(i),
-      this->neighbor_index(i),
-      this->dof_handler);
-
-#ifdef DEBUG
-  if (q.state() != IteratorState::past_the_end)
-    Assert(q->used(), ExcInternalError());
-#endif
-  return q;
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline TriaIterator<
-  DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::child(
-  const unsigned int i) const
-{
-  TriaIterator<DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
-    q(this->tria, this->level() + 1, this->child_index(i), this->dof_handler);
-
-#ifdef DEBUG
-  if (q.state() != IteratorState::past_the_end)
-    Assert(q->used(), ExcInternalError());
-#endif
-  return q;
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline boost::container::small_vector<
-  TriaIterator<DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>,
-  GeometryInfo<dimension_>::max_children_per_cell>
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  child_iterators() const
-{
-  boost::container::small_vector<
-    TriaIterator<
-      DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>,
-    GeometryInfo<dimension_>::max_children_per_cell>
-    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 dimension_, int space_dimension_, bool level_dof_access>
-inline TriaIterator<
-  DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::parent() const
-{
-  TriaIterator<DoFCellAccessor<dimension_, space_dimension_, level_dof_access>>
-    q(this->tria, this->level() - 1, this->parent_index(), this->dof_handler);
-
-  return q;
-}
-
-
-
-namespace internal
-{
-  namespace DoFCellAccessorImplementation
-  {
-    template <int dim, int spacedim, bool level_dof_access>
-    inline TriaIterator<
-      dealii::DoFAccessor<dim - 1, dim, spacedim, level_dof_access>>
-    get_face(
-      const dealii::DoFCellAccessor<dim, spacedim, level_dof_access> &cell,
-      const unsigned int                                              i,
-      const std::integral_constant<int, 1>)
-    {
-      dealii::DoFAccessor<0, dim, spacedim, level_dof_access> a(
-        &cell.get_triangulation(),
-        ((i == 0) && cell.at_boundary(0) ?
-           dealii::TriaAccessor<0, 1, spacedim>::left_vertex :
-           ((i == 1) && cell.at_boundary(1) ?
-              dealii::TriaAccessor<0, 1, spacedim>::right_vertex :
-              dealii::TriaAccessor<0, 1, spacedim>::interior_vertex)),
-        cell.vertex_index(i),
-        &cell.get_dof_handler());
-      return dealii::TriaIterator<
-        dealii::DoFAccessor<0, dim, spacedim, level_dof_access>>(a);
-    }
-
-
-    template <int dim, int spacedim, bool level_dof_access>
-    inline TriaIterator<
-      dealii::DoFAccessor<dim - 1, dim, spacedim, level_dof_access>>
-    get_face(
-      const dealii::DoFCellAccessor<dim, spacedim, level_dof_access> &cell,
-      const unsigned int                                              i,
-      const std::integral_constant<int, 2>)
-    {
-      return cell.line(i);
-    }
-
-
-    template <int dim, int spacedim, bool level_dof_access>
-    inline TriaIterator<
-      dealii::DoFAccessor<dim - 1, dim, spacedim, level_dof_access>>
-    get_face(
-      const dealii::DoFCellAccessor<dim, spacedim, level_dof_access> &cell,
-      const unsigned int                                              i,
-      const std::integral_constant<int, 3>)
-    {
-      return cell.quad(i);
-    }
-  } // namespace DoFCellAccessorImplementation
-} // namespace internal
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline typename DoFCellAccessor<dimension_,
-                                space_dimension_,
-                                level_dof_access>::face_iterator
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::face(
-  const unsigned int i) const
-{
-  AssertIndexRange(i, this->n_faces());
-
-  return dealii::internal::DoFCellAccessorImplementation::get_face(
-    *this, i, std::integral_constant<int, dimension_>());
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline boost::container::small_vector<
-  typename DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-    face_iterator,
-  GeometryInfo<dimension_>::faces_per_cell>
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  face_iterators() const
-{
-  boost::container::small_vector<
-    typename DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-      face_iterator,
-    GeometryInfo<dimension_>::faces_per_cell>
-    face_iterators(this->n_faces());
-
-  for (const unsigned int i : this->face_indices())
-    face_iterators[i] =
-      dealii::internal::DoFCellAccessorImplementation::get_face(
-        *this, i, std::integral_constant<int, dimension_>());
-
-  return face_iterators;
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  get_active_or_mg_dof_indices(
-    std::vector<types::global_dof_index> &dof_indices) const
-{
-  if (level_dof_access)
-    get_mg_dof_indices(dof_indices);
-  else
-    get_dof_indices(dof_indices);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <class InputVector, typename number>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::get_dof_values(
-  const InputVector &values,
-  Vector<number>    &local_values) const
-{
-  get_dof_values(values, local_values.begin(), local_values.end());
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <typename Number, typename ForwardIterator>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::get_dof_values(
-  const ReadVector<Number> &values,
-  ForwardIterator           local_values_begin,
-  ForwardIterator           local_values_end) const
-{
-  (void)local_values_end;
-  Assert(this->is_artificial() == false,
-         ExcMessage("Can't ask for DoF indices on artificial cells."));
-  Assert(this->is_active(), ExcMessage("Cell must be active."));
-  Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
-
-  Assert(static_cast<unsigned int>(local_values_end - local_values_begin) ==
-           this->get_fe().n_dofs_per_cell(),
-         typename DoFCellAccessor::ExcVectorDoesNotMatch());
-  Assert(values.size() == this->get_dof_handler().n_dofs(),
-         typename DoFCellAccessor::ExcVectorDoesNotMatch());
-
-  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
-    dof_indices(this->get_fe().n_dofs_per_cell());
-  internal::DoFAccessorImplementation::get_cell_dof_indices(
-    *this, dof_indices, this->active_fe_index());
-
-  boost::container::small_vector<Number, 27> values_temp(local_values_end -
-                                                         local_values_begin);
-  auto view = make_array_view(values_temp.begin(), values_temp.end());
-  values.extract_subvector_to(make_array_view(dof_indices.begin(),
-                                              dof_indices.end()),
-                              view);
-  using view_type = std::remove_reference_t<decltype(*local_values_begin)>;
-  ArrayView<view_type> values_view2(&*local_values_begin,
-                                    local_values_end - local_values_begin);
-  std::copy(values_temp.begin(), values_temp.end(), values_view2.begin());
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <class InputVector, typename ForwardIterator>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::get_dof_values(
-  const AffineConstraints<typename InputVector::value_type> &constraints,
-  const InputVector                                         &values,
-  ForwardIterator                                            local_values_begin,
-  ForwardIterator local_values_end) const
-{
-  Assert(this->is_artificial() == false,
-         ExcMessage("Can't ask for DoF indices on artificial cells."));
-  Assert(this->is_active(), ExcMessage("Cell must be active."));
-
-  Assert(static_cast<unsigned int>(local_values_end - local_values_begin) ==
-           this->get_fe().n_dofs_per_cell(),
-         typename DoFCellAccessor::ExcVectorDoesNotMatch());
-  Assert(values.size() == this->get_dof_handler().n_dofs(),
-         typename DoFCellAccessor::ExcVectorDoesNotMatch());
-
-
-  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
-    dof_indices(this->get_fe().n_dofs_per_cell());
-  internal::DoFAccessorImplementation::get_cell_dof_indices(
-    *this, dof_indices, this->active_fe_index());
-
-  constraints.get_dof_values(values,
-                             dof_indices.data(),
-                             local_values_begin,
-                             local_values_end);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <class OutputVector, typename number>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::set_dof_values(
-  const Vector<number> &local_values,
-  OutputVector         &values) const
-{
-  Assert(this->is_artificial() == false,
-         ExcMessage("Can't ask for DoF indices on artificial cells."));
-  Assert(this->is_active(), ExcMessage("Cell must be active."));
-
-  Assert(static_cast<unsigned int>(local_values.size()) ==
-           this->get_fe().n_dofs_per_cell(),
-         typename DoFCellAccessor::ExcVectorDoesNotMatch());
-  Assert(values.size() == this->get_dof_handler().n_dofs(),
-         typename DoFCellAccessor::ExcVectorDoesNotMatch());
-
-
-  Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
-  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
-    dof_indices(this->get_fe().n_dofs_per_cell());
-  internal::DoFAccessorImplementation::get_cell_dof_indices(
-    *this, dof_indices, this->active_fe_index());
-
-  for (unsigned int i = 0; i < this->get_fe().n_dofs_per_cell(); ++i)
-    internal::ElementAccess<OutputVector>::set(local_values(i),
-                                               dof_indices[i],
-                                               values);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline const FiniteElement<dimension_, space_dimension_> &
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::get_fe() const
-{
-  Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           this->is_active(),
-         ExcMessage(
-           "For DoFHandler objects in hp-mode, finite elements are only "
-           "associated with active cells. Consequently, you can not ask "
-           "for the active finite element on cells with children."));
-
-  const auto &fe = this->dof_handler->get_fe(active_fe_index());
-
-  Assert(this->reference_cell() == fe.reference_cell(),
-         internal::ExcNonMatchingReferenceCellTypes(this->reference_cell(),
-                                                    fe.reference_cell()));
-
-  return fe;
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline types::fe_index
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  active_fe_index() const
-{
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           this->is_active(),
-         ExcMessage(
-           "You can not ask for the active FE index on a cell that has "
-           "children because no degrees of freedom are assigned "
-           "to this cell and, consequently, no finite element "
-           "is associated with it."));
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           (this->is_locally_owned() || this->is_ghost()),
-         ExcMessage("You can only query active FE index information on cells "
-                    "that are either locally owned or (after distributing "
-                    "degrees of freedom) are ghost cells."));
-
-  return dealii::internal::DoFCellAccessorImplementation::Implementation::
-    active_fe_index(*this);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  set_active_fe_index(const types::fe_index i) const
-{
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           this->is_active(),
-         ExcMessage("You can not set the active FE index on a cell that has "
-                    "children because no degrees of freedom will be assigned "
-                    "to this cell."));
-
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           this->is_locally_owned(),
-         ExcMessage("You can only set active FE index information on cells "
-                    "that are locally owned. On ghost cells, this information "
-                    "will automatically be propagated from the owning process "
-                    "of that cell, and there is no information at all on "
-                    "artificial cells."));
-
-  dealii::internal::DoFCellAccessorImplementation::Implementation::
-    set_active_fe_index(*this, i);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline const FiniteElement<dimension_, space_dimension_> &
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::get_future_fe()
-  const
-{
-  Assert(this->dof_handler != nullptr, typename BaseClass::ExcInvalidObject());
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           this->is_active(),
-         ExcMessage(
-           "For DoFHandler objects in hp-mode, finite elements are only "
-           "associated with active cells. Consequently, you can not ask "
-           "for the future finite element on cells with children."));
-
-  return this->dof_handler->get_fe(future_fe_index());
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline types::fe_index
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  future_fe_index() const
-{
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           (this->has_children() == false),
-         ExcMessage(
-           "You can not ask for the future FE index on a cell that has "
-           "children because no degrees of freedom are assigned "
-           "to this cell and, consequently, no finite element "
-           "is associated with it."));
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           (this->is_locally_owned()),
-         ExcMessage("You can only query future FE index information on cells "
-                    "that are locally owned."));
-
-  return dealii::internal::DoFCellAccessorImplementation::Implementation::
-    future_fe_index(*this);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  set_future_fe_index(const types::fe_index i) const
-{
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           (this->has_children() == false),
-         ExcMessage("You can not set the future FE index on a cell that has "
-                    "children because no degrees of freedom will be assigned "
-                    "to this cell."));
-
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           this->is_locally_owned(),
-         ExcMessage("You can only set future FE index information on cells "
-                    "that are locally owned."));
-
-  dealii::internal::DoFCellAccessorImplementation::Implementation::
-    set_future_fe_index(*this, i);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline bool
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  future_fe_index_set() const
-{
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           (this->has_children() == false),
-         ExcMessage(
-           "You can not ask for the future FE index on a cell that has "
-           "children because no degrees of freedom are assigned "
-           "to this cell and, consequently, no finite element "
-           "is associated with it."));
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           (this->is_locally_owned()),
-         ExcMessage("You can only query future FE index information on cells "
-                    "that are locally owned."));
-
-  return dealii::internal::DoFCellAccessorImplementation::Implementation::
-    future_fe_index_set(*this);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  clear_future_fe_index() const
-{
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           (this->has_children() == false),
-         ExcMessage(
-           "You can not ask for the future FE index on a cell that has "
-           "children because no degrees of freedom are assigned "
-           "to this cell and, consequently, no finite element "
-           "is associated with it."));
-  Assert((this->dof_handler->hp_capability_enabled == false) ||
-           (this->is_locally_owned()),
-         ExcMessage("You can only query future FE index information on cells "
-                    "that are locally owned."));
-
-  dealii::internal::DoFCellAccessorImplementation::Implementation::
-    clear_future_fe_index(*this);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <typename number, typename OutputVector>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  distribute_local_to_global(const Vector<number> &local_source,
-                             OutputVector         &global_destination) const
-{
-  this->distribute_local_to_global(local_source.begin(),
-                                   local_source.end(),
-                                   global_destination);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <typename ForwardIterator, typename OutputVector>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  distribute_local_to_global(ForwardIterator local_source_begin,
-                             ForwardIterator local_source_end,
-                             OutputVector   &global_destination) const
-{
-  Assert(this->dof_handler != nullptr,
-         (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
-  Assert(static_cast<unsigned int>(local_source_end - local_source_begin) ==
-           this->get_fe().n_dofs_per_cell(),
-         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
-  Assert(this->dof_handler->n_dofs() == global_destination.size(),
-         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
-
-  Assert(!this->has_children(), ExcMessage("Cell must be active"));
-
-  Assert(
-    internal::ArrayViewHelper::is_contiguous(local_source_begin,
-                                             local_source_end),
-    ExcMessage(
-      "This function can not be called with iterator types that do not point to contiguous memory."));
-
-  const unsigned int n_dofs = local_source_end - local_source_begin;
-
-  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
-    dof_indices(n_dofs);
-  internal::DoFAccessorImplementation::get_cell_dof_indices(
-    *this, dof_indices, this->active_fe_index());
-
-  // distribute cell vector
-  global_destination.add(n_dofs, dof_indices.data(), &(*local_source_begin));
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <typename ForwardIterator, typename OutputVector>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  distribute_local_to_global(
-    const AffineConstraints<typename OutputVector::value_type> &constraints,
-    ForwardIterator local_source_begin,
-    ForwardIterator local_source_end,
-    OutputVector   &global_destination) const
-{
-  Assert(this->dof_handler != nullptr,
-         (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
-  Assert(local_source_end - local_source_begin ==
-           this->get_fe().n_dofs_per_cell(),
-         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
-  Assert(this->dof_handler->n_dofs() == global_destination.size(),
-         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
-
-  Assert(!this->has_children(), ExcMessage("Cell must be active."));
-
-  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
-    dof_indices(this->get_fe().n_dofs_per_cell());
-  internal::DoFAccessorImplementation::get_cell_dof_indices(
-    *this, dof_indices, this->active_fe_index());
-
-  // distribute cell vector
-  constraints.distribute_local_to_global(local_source_begin,
-                                         local_source_end,
-                                         dof_indices.data(),
-                                         global_destination);
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <typename number, typename OutputMatrix>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  distribute_local_to_global(const FullMatrix<number> &local_source,
-                             OutputMatrix             &global_destination) const
-{
-  Assert(this->dof_handler != nullptr,
-         (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
-  Assert(local_source.m() == this->get_fe().n_dofs_per_cell(),
-         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
-  Assert(local_source.n() == this->get_fe().n_dofs_per_cell(),
-         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
-  Assert(this->dof_handler->n_dofs() == global_destination.m(),
-         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
-  Assert(this->dof_handler->n_dofs() == global_destination.n(),
-         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
-
-  Assert(!this->has_children(), ExcMessage("Cell must be active."));
-
-  const unsigned int n_dofs = local_source.m();
-
-  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
-    dof_indices(n_dofs);
-  internal::DoFAccessorImplementation::get_cell_dof_indices(
-    *this, dof_indices, this->active_fe_index());
-
-  // distribute cell matrix
-  for (unsigned int i = 0; i < n_dofs; ++i)
-    global_destination.add(dof_indices[i],
-                           n_dofs,
-                           dof_indices.data(),
-                           &local_source(i, 0));
-}
-
-
-
-template <int dimension_, int space_dimension_, bool level_dof_access>
-template <typename number, typename OutputMatrix, typename OutputVector>
-inline void
-DoFCellAccessor<dimension_, space_dimension_, level_dof_access>::
-  distribute_local_to_global(const FullMatrix<number> &local_matrix,
-                             const Vector<number>     &local_vector,
-                             OutputMatrix             &global_matrix,
-                             OutputVector             &global_vector) const
-{
-  Assert(this->dof_handler != nullptr,
-         (typename std::decay_t<decltype(*this)>::ExcInvalidObject()));
-  Assert(local_matrix.m() == this->get_fe().n_dofs_per_cell(),
-         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
-  Assert(local_matrix.n() == this->get_fe().n_dofs_per_cell(),
-         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
-  Assert(this->dof_handler->n_dofs() == global_matrix.m(),
-         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
-  Assert(this->dof_handler->n_dofs() == global_matrix.n(),
-         (typename std::decay_t<decltype(*this)>::ExcMatrixDoesNotMatch()));
-  Assert(local_vector.size() == this->get_fe().n_dofs_per_cell(),
-         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
-  Assert(this->dof_handler->n_dofs() == global_vector.size(),
-         (typename std::decay_t<decltype(*this)>::ExcVectorDoesNotMatch()));
-
-  Assert(!this->has_children(), ExcMessage("Cell must be active."));
-
-  const unsigned int n_dofs = this->get_fe().n_dofs_per_cell();
-  internal::DoFAccessorImplementation::Implementation::dof_index_vector_type
-    dof_indices(n_dofs);
-  internal::DoFAccessorImplementation::get_cell_dof_indices(
-    *this, dof_indices, this->active_fe_index());
-
-  // distribute cell matrices
-  for (unsigned int i = 0; i < n_dofs; ++i)
-    {
-      global_matrix.add(dof_indices[i],
-                        n_dofs,
-                        dof_indices.data(),
-                        &local_matrix(i, 0));
-      global_vector(dof_indices[i]) += local_vector(i);
-    }
-}
-
-
-
-DEAL_II_NAMESPACE_CLOSE
+DEAL_II_WARNING(
+  "The use of this header file is deprecated. Just include <deal.II/dofs/dof_accessor.h>.")
 
 #endif

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.