fe_evaluation_data.h

Go to the documentation of this file.

// -----------------------------------------------------------------------------
//
// SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception OR LGPL-2.1-or-later
// Copyright (C) 2021 - 2026 by the deal.II authors
//
// This file is part of the deal.II library.
//
// Detailed license information governing the source code and contributions
// can be found in LICENSE.md and CONTRIBUTING.md at the top level directory.
//
// -----------------------------------------------------------------------------
 
 
#ifndef dealii_matrix_free_fe_evaluation_data_h
#define dealii_matrix_free_fe_evaluation_data_h
 
 
#include <deal.II/base/config.h>
 
#include <deal.II/base/array_view.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/geometry_info.h>
#include <deal.II/base/observer_pointer.h>
#include <deal.II/base/signaling_nan.h>
#include <deal.II/base/std_cxx20/iota_view.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/template_constraints.h>
#include <deal.II/base/tensor.h>
#include <deal.II/base/vectorization.h>
 
#include <deal.II/grid/reference_cell.h>
 
#include <deal.II/matrix_free/dof_info.h>
#include <deal.II/matrix_free/mapping_info_storage.h>
#include <deal.II/matrix_free/shape_info.h>
 
 
DEAL_II_NAMESPACE_OPEN
 
 
 
namespace internal
{
  DeclException0(ExcAccessToUninitializedField);
 
  DeclException1(
    ExcMatrixFreeAccessToUninitializedMappingField,
    std::string,
    << "You are requesting information from an FEEvaluation/FEFaceEvaluation "
    << "object for which this kind of information has not been computed. What "
    << "information these objects compute is determined by the update_* flags "
    << "you pass to MatrixFree::reinit() via MatrixFree::AdditionalData. "
    << "Here, the operation you are attempting requires the <" << arg1
    << "> flag to be set, but it was apparently not specified "
    << "upon initialization.");
} // namespace internal
 
// forward declarations
template <int dim,
          int fe_degree,
          int n_q_points_1d            = fe_degree + 1,
          int n_components_            = 1,
          typename Number              = double,
          typename VectorizedArrayType = VectorizedArray<Number>>
class FEEvaluation;
 
template <int dim,
          int n_components_,
          typename Number,
          bool     = false,
          typename = VectorizedArray<Number>>
class FEEvaluationBase;
 
namespace internal
{
  namespace MatrixFreeFunctions
  {
    template <int, typename>
    class MappingDataOnTheFly;
  }
} // namespace internal
 
 
template <int dim, typename Number, bool is_face>
class FEEvaluationData
{
  using ShapeInfoType = internal::MatrixFreeFunctions::ShapeInfo<
    typename internal::VectorizedArrayTrait<Number>::value_type>;
  using MappingInfoStorageType = internal::MatrixFreeFunctions::
    MappingInfoStorage<(is_face ? dim - 1 : dim), dim, Number>;
  using DoFInfo = internal::MatrixFreeFunctions::DoFInfo;
 
public:
  static constexpr unsigned int dimension = dim;
 
  using NumberType = Number;
  using ScalarNumber =
    typename internal::VectorizedArrayTrait<Number>::value_type;
  using shape_info_number_type = ScalarNumber;
 
  static constexpr unsigned int n_lanes =
    internal::VectorizedArrayTrait<Number>::width();
 
  FEEvaluationData(const ShapeInfoType &shape_info,
                   const bool           is_interior_face = true);
 
  FEEvaluationData(const FEEvaluationData &other) = default;
 
  FEEvaluationData &
  operator=(const FEEvaluationData &other);
 
  virtual ~FEEvaluationData() = default;
 
  void
  set_data_pointers(AlignedVector<Number> *scratch_data,
                    const unsigned int     n_components);
 
  void
  reinit_face(
    const internal::MatrixFreeFunctions::FaceToCellTopology<n_lanes> &face);
 
  ArrayView<Number>
  get_scratch_data() const;
 
  Number
  JxW(const unsigned int q_point) const;
 
  Point<dim, Number>
  quadrature_point(const unsigned int q) const;
 
  Tensor<2, dim, Number>
  inverse_jacobian(const unsigned int q_point) const;
 
  Tensor<1, dim, Number>
  normal_vector(const unsigned int q_point) const;
 
  DEAL_II_DEPRECATED_WITH_COMMENT("Use normal_vector() instead.")
  Tensor<1, dim, Number>
  get_normal_vector(const unsigned int q_point) const;
 
  const Number *
  begin_dof_values() const;
 
  Number *
  begin_dof_values();
 
  const Number *
  begin_values() const;
 
  Number *
  begin_values();
 
  const Number *
  begin_gradients() const;
 
  Number *
  begin_gradients();
 
  const Number *
  begin_hessians() const;
 
  Number *
  begin_hessians();
 
  unsigned int
  get_mapping_data_index_offset() const;
 
  internal::MatrixFreeFunctions::GeometryType
  get_cell_type() const;
 
  const ShapeInfoType &
  get_shape_info() const;
 
  const internal::MatrixFreeFunctions::DoFInfo &
  get_dof_info() const;
 
  const std::vector<unsigned int> &
  get_internal_dof_numbering() const;
 
  unsigned int
  get_quadrature_index() const;
 
  unsigned int
  get_current_cell_index() const;
 
  unsigned int
  get_active_fe_index() const;
 
  unsigned int
  get_active_quadrature_index() const;
 
  unsigned int
  get_first_selected_component() const;
 
  std::uint8_t
  get_face_no(const unsigned int v = 0) const;
 
  unsigned int
  get_subface_index() const;
 
  std::uint8_t
  get_face_orientation(const unsigned int v = 0) const;
 
  internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex
  get_dof_access_index() const;
 
  bool
  is_interior_face() const;
 
  const std::array<unsigned int, n_lanes> &
  get_cell_ids() const
  {
    // implemented inline to avoid compilation problems on Windows
    if constexpr (running_in_debug_mode())
      {
        Assert(is_reinitialized, ExcNotInitialized());
      }
    return cell_ids;
  }
 
  const std::array<unsigned int, n_lanes> &
  get_face_ids() const
  {
    // implemented inline to avoid compilation problems on Windows
    if constexpr (running_in_debug_mode())
      {
        Assert(is_reinitialized && is_face, ExcNotInitialized());
      }
    return face_ids;
  }
 
  unsigned int
  get_cell_or_face_batch_id() const
  {
    // implemented inline to avoid compilation problems on Windows
    if constexpr (running_in_debug_mode())
      {
        Assert(is_reinitialized, ExcNotInitialized());
      }
 
    return cell;
  }
 
  const std::array<unsigned int, n_lanes> &
  get_cell_or_face_ids() const
  {
    // implemented inline to avoid compilation problems on Windows
    if constexpr (running_in_debug_mode())
      {
        Assert(is_reinitialized, ExcNotInitialized());
      }
 
    if (!is_face || dof_access_index ==
                      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
      return cell_ids;
    else
      return face_ids;
  }
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  quadrature_point_indices() const;
 
  template <typename T>
  T
  read_cell_data(const AlignedVector<T> &array) const;
 
  template <typename T, int D>
  void
  read_cell_data(const Table<D, T> &src, Table<D - 1, T> &dst) const;
 
  template <typename T>
  void
  set_cell_data(AlignedVector<T> &array, const T &value) const;
 
  template <typename T>
  T
  read_face_data(const AlignedVector<T> &array) const;
 
  template <typename T>
  void
  set_face_data(AlignedVector<T> &array, const T &value) const;
 
  struct InitializationData
  {
    const ShapeInfoType          *shape_info;
    const DoFInfo                *dof_info;
    const MappingInfoStorageType *mapping_data;
    unsigned int                  active_fe_index;
    unsigned int                  active_quad_index;
    const typename MappingInfoStorageType::QuadratureDescriptor *descriptor;
  };
 
protected:
  FEEvaluationData(const InitializationData &initialization_data,
                   const bool                is_interior_face,
                   const unsigned int        quadrature_index,
                   const unsigned int        first_selected_component);
 
  FEEvaluationData(
    const std::shared_ptr<
      internal::MatrixFreeFunctions::MappingDataOnTheFly<dim, Number>>
                      &mapping_data,
    const unsigned int n_fe_components,
    const unsigned int first_selected_component);
 
  const ShapeInfoType *data;
 
  const DoFInfo *dof_info;
 
  const MappingInfoStorageType *mapping_data;
 
  const unsigned int quadrature_index;
 
  const unsigned int n_fe_components;
 
  const unsigned int first_selected_component;
 
  const unsigned int active_fe_index;
 
  const unsigned int active_quad_index;
 
  const typename MappingInfoStorageType::QuadratureDescriptor *descriptor;
 
  const unsigned int n_quadrature_points;
 
  const Point<dim, Number> *quadrature_points;
 
  const Tensor<2, dim, Number> *jacobian;
 
  const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, Number>>
    *jacobian_gradients;
 
  const Tensor<1, dim *(dim + 1) / 2, Tensor<1, dim, Number>>
    *jacobian_gradients_non_inverse;
 
  const Number *J_value;
 
  const Tensor<1, dim, Number> *normal_vectors;
 
  const Tensor<1, dim, Number> *normal_x_jacobian;
 
  const ScalarNumber *quadrature_weights;
 
  mutable ArrayView<Number> scratch_data;
 
  Number *values_dofs;
 
  Number *values_quad;
 
  Number *values_from_gradients_quad;
 
  Number *gradients_quad;
 
  Number *gradients_from_hessians_quad;
 
  Number *hessians_quad;
 
  bool is_reinitialized;
 
  bool dof_values_initialized;
 
  bool values_quad_initialized;
 
  bool gradients_quad_initialized;
 
  bool hessians_quad_initialized;
 
  bool values_quad_submitted;
 
  bool gradients_quad_submitted;
 
  bool hessians_quad_submitted;
 
  unsigned int cell;
 
  bool interior_face;
 
  internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex dof_access_index;
 
  std::array<std::uint8_t, n_lanes> face_numbers;
 
  std::array<types::geometric_orientation, n_lanes> face_orientations;
 
  unsigned int subface_index;
 
  internal::MatrixFreeFunctions::GeometryType cell_type;
 
  std::array<unsigned int, n_lanes> cell_ids;
 
  std::array<unsigned int, n_lanes> face_ids;
 
  std::shared_ptr<
    internal::MatrixFreeFunctions::MappingDataOnTheFly<dim, Number>>
    mapped_geometry;
 
  bool divergence_is_requested;
 
  // Make FEEvaluation and FEEvaluationBase objects friends for access to
  // protected member mapped_geometry.
  template <int, int, typename, bool, typename>
  friend class FEEvaluationBase;
 
  template <int, int, int, int, typename, typename>
  friend class FEEvaluation;
};
 
 
template <int dim,
          typename Number,
          bool is_face,
          typename VectorizedArrayType = VectorizedArray<Number>>
using FEEvaluationBaseData =
  FEEvaluationData<dim, VectorizedArrayType, is_face>;
 
 
/*----------------------- Inline functions ----------------------------------*/
 
#ifndef DOXYGEN
 
template <int dim, typename Number, bool is_face>
inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
  const ShapeInfoType &shape_info,
  const bool           is_interior_face)
  : FEEvaluationData(
      InitializationData{&shape_info, nullptr, nullptr, 0, 0, nullptr},
      is_interior_face,
      0,
      0)
{}
 
 
template <int dim, typename Number, bool is_face>
inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
  const InitializationData &initialization_data,
  const bool                is_interior_face,
  const unsigned int        quadrature_index,
  const unsigned int        first_selected_component)
  : data(initialization_data.shape_info)
  , dof_info(initialization_data.dof_info)
  , mapping_data(initialization_data.mapping_data)
  , quadrature_index(quadrature_index)
  , n_fe_components(dof_info != nullptr ? dof_info->start_components.back() : 0)
  , first_selected_component(first_selected_component)
  , active_fe_index(initialization_data.active_fe_index)
  , active_quad_index(initialization_data.active_quad_index)
  , descriptor(initialization_data.descriptor)
  , n_quadrature_points(descriptor == nullptr ?
                          (is_face ? data->n_q_points_face : data->n_q_points) :
                          descriptor->n_q_points)
  , quadrature_points(nullptr)
  , jacobian(nullptr)
  , jacobian_gradients(nullptr)
  , jacobian_gradients_non_inverse(nullptr)
  , J_value(nullptr)
  , normal_vectors(nullptr)
  , normal_x_jacobian(nullptr)
  , quadrature_weights(
      descriptor != nullptr ? descriptor->quadrature_weights.begin() : nullptr)
#  ifdef DEBUG
  , is_reinitialized(false)
  , dof_values_initialized(false)
  , values_quad_initialized(false)
  , gradients_quad_initialized(false)
  , hessians_quad_initialized(false)
  , values_quad_submitted(false)
  , gradients_quad_submitted(false)
  , hessians_quad_submitted(false)
#  endif
  , cell(numbers::invalid_unsigned_int)
  , interior_face(is_interior_face)
  , dof_access_index(
      is_face ?
        (is_interior_face ?
           internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
           internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
        internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
  , subface_index(0)
  , cell_type(internal::MatrixFreeFunctions::general)
  , divergence_is_requested(false)
{
  Assert(!data->data.empty(), ExcInternalError());
}
 
 
 
template <int dim, typename Number, bool is_face>
inline FEEvaluationData<dim, Number, is_face>::FEEvaluationData(
  const std::shared_ptr<
    internal::MatrixFreeFunctions::MappingDataOnTheFly<dim, Number>>
                    &mapped_geometry,
  const unsigned int n_fe_components,
  const unsigned int first_selected_component)
  : data(nullptr)
  , dof_info(nullptr)
  , mapping_data(&mapped_geometry->get_data_storage())
  , quadrature_index(numbers::invalid_unsigned_int)
  , n_fe_components(n_fe_components)
  , first_selected_component(first_selected_component)
  , active_fe_index(numbers::invalid_unsigned_int)
  , active_quad_index(numbers::invalid_unsigned_int)
  , descriptor(&mapped_geometry->get_data_storage().descriptor[0])
  , n_quadrature_points(descriptor->n_q_points)
  , quadrature_points(
      mapped_geometry->get_data_storage().quadrature_points.begin())
  , jacobian(mapped_geometry->get_data_storage().jacobians[0].begin())
  , jacobian_gradients(
      mapped_geometry->get_data_storage().jacobian_gradients[0].begin())
  , jacobian_gradients_non_inverse(mapped_geometry->get_data_storage()
                                     .jacobian_gradients_non_inverse[0]
                                     .begin())
  , J_value(mapped_geometry->get_data_storage().JxW_values.begin())
  , normal_vectors(nullptr)
  , normal_x_jacobian(nullptr)
  , quadrature_weights(descriptor->quadrature_weights.begin())
  , cell(0)
  , cell_type(internal::MatrixFreeFunctions::general)
  , interior_face(true)
  , dof_access_index(internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
  , mapped_geometry(mapped_geometry)
  , is_reinitialized(false)
  , divergence_is_requested(false)
{}
 
 
 
template <int dim, typename Number, bool is_face>
inline FEEvaluationData<dim, Number, is_face> &
FEEvaluationData<dim, Number, is_face>::operator=(const FEEvaluationData &other)
{
  AssertDimension(quadrature_index, other.quadrature_index);
  AssertDimension(n_fe_components, other.n_fe_components);
  AssertDimension(first_selected_component, other.first_selected_component);
  AssertDimension(active_fe_index, other.active_fe_index);
  AssertDimension(active_quad_index, other.active_quad_index);
  AssertDimension(n_quadrature_points, descriptor->n_q_points);
 
  data                           = other.data;
  dof_info                       = other.dof_info;
  mapping_data                   = other.mapping_data;
  descriptor                     = other.descriptor;
  jacobian                       = nullptr;
  J_value                        = nullptr;
  normal_vectors                 = nullptr;
  normal_x_jacobian              = nullptr;
  jacobian_gradients             = nullptr;
  jacobian_gradients_non_inverse = nullptr;
  quadrature_points              = nullptr;
  quadrature_weights             = other.quadrature_weights;
 
  if constexpr (running_in_debug_mode())
    {
      is_reinitialized           = false;
      dof_values_initialized     = false;
      values_quad_initialized    = false;
      gradients_quad_initialized = false;
      hessians_quad_initialized  = false;
      hessians_quad_submitted    = false;
      values_quad_submitted      = false;
      gradients_quad_submitted   = false;
    }
 
  cell          = numbers::invalid_unsigned_int;
  interior_face = other.is_interior_face();
  dof_access_index =
    is_face ?
      (is_interior_face() ?
         internal::MatrixFreeFunctions::DoFInfo::dof_access_face_interior :
         internal::MatrixFreeFunctions::DoFInfo::dof_access_face_exterior) :
      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell;
  face_numbers[0]         = 0;
  face_orientations[0]    = numbers::default_geometric_orientation;
  subface_index           = 0;
  cell_type               = internal::MatrixFreeFunctions::general;
  divergence_is_requested = false;
 
  return *this;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline void
FEEvaluationData<dim, Number, is_face>::set_data_pointers(
  AlignedVector<Number> *scratch_data_array,
  const unsigned int     n_components)
{
  Assert(scratch_data_array != nullptr, ExcInternalError());
  Assert(data != nullptr, ExcInternalError());
  Assert(!data->data.empty(), ExcInternalError());
 
 
  const unsigned int tensor_dofs_per_component =
    Utilities::fixed_power<dim>(data->data.front().fe_degree + 1);
  const unsigned int dofs_per_component = data->dofs_per_component_on_cell;
 
  const unsigned int size_scratch_data =
    std::max(tensor_dofs_per_component + 1, dofs_per_component) * n_components *
      4 +
    2 * n_quadrature_points;
  const unsigned int size_data_arrays =
    n_components * dofs_per_component +
    (n_components * ((dim * (dim + 1)) / 2 + 2 * dim + 2) *
     n_quadrature_points);
 
  // include 12 extra fields to insert some padding between values, gradients
  // and hessians, which helps to reduce the probability of cache conflicts
  const unsigned int allocated_size = size_scratch_data + size_data_arrays + 12;
  if constexpr (running_in_debug_mode())
    {
      scratch_data_array->resize(allocated_size);
      scratch_data_array->fill(Number(numbers::signaling_nan<ScalarNumber>()));
    }
  else
    {
      scratch_data_array->resize_fast(allocated_size);
    }
  scratch_data.reinit(scratch_data_array->begin() + size_data_arrays + 12,
                      size_scratch_data);
 
  // set the pointers to the correct position in the data array
  values_dofs = scratch_data_array->begin();
  values_quad =
    scratch_data_array->begin() + 4 + n_components * dofs_per_component;
  values_from_gradients_quad =
    scratch_data_array->begin() + 6 +
    n_components * (dofs_per_component + n_quadrature_points);
  gradients_quad =
    scratch_data_array->begin() + 8 +
    n_components * (dofs_per_component + 2 * n_quadrature_points);
  gradients_from_hessians_quad =
    scratch_data_array->begin() + 12 +
    n_components * (dofs_per_component + (dim + 2) * n_quadrature_points);
  hessians_quad =
    scratch_data_array->begin() + 12 +
    n_components * (dofs_per_component + (2 * dim + 2) * n_quadrature_points);
}
 
 
 
template <int dim, typename Number, bool is_face>
inline void
FEEvaluationData<dim, Number, is_face>::reinit_face(
  const internal::MatrixFreeFunctions::FaceToCellTopology<n_lanes> &face)
{
  Assert(is_face == true,
         ExcMessage("Faces can only be set if the is_face template parameter "
                    "is true"));
  face_numbers[0] =
    (is_interior_face() ? face.interior_face_no : face.exterior_face_no);
  subface_index = is_interior_face() == true ?
                    GeometryInfo<dim>::max_children_per_cell :
                    face.subface_index;
 
  // First check if interior or exterior cell has non-standard orientation
  // (i.e. the third bit is one or not). Then set zero if this cell has
  // standard-orientation else copy the first three bits
  // (which is equivalent to modulo 8). See also the documentation of
  // internal::MatrixFreeFunctions::FaceToCellTopology::face_orientation.
  face_orientations[0] = (is_interior_face() == (face.face_orientation >= 8)) ?
                           (face.face_orientation % 8) :
                           numbers::default_geometric_orientation;
 
  if (is_interior_face())
    cell_ids = face.cells_interior;
  else
    cell_ids = face.cells_exterior;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
FEEvaluationData<dim, Number, is_face>::normal_vector(
  const unsigned int q_point) const
{
  AssertIndexRange(q_point, n_quadrature_points);
  Assert(normal_vectors != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_normal_vectors"));
  if (cell_type <= internal::MatrixFreeFunctions::flat_faces)
    return normal_vectors[0];
  else
    return normal_vectors[q_point];
}
 
 
 
// This function is deprecated.
template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Tensor<1, dim, Number>
FEEvaluationData<dim, Number, is_face>::get_normal_vector(
  const unsigned int q_point) const
{
  return normal_vector(q_point);
}
 
 
 
template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Number
FEEvaluationData<dim, Number, is_face>::JxW(const unsigned int q_point) const
{
  AssertIndexRange(q_point, n_quadrature_points);
  Assert(J_value != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_values|update_gradients"));
  if (cell_type <= internal::MatrixFreeFunctions::affine)
    {
      Assert(quadrature_weights != nullptr, ExcNotInitialized());
      return J_value[0] * quadrature_weights[q_point];
    }
  else
    return J_value[q_point];
}
 
 
 
template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Point<dim, Number>
FEEvaluationData<dim, Number, is_face>::quadrature_point(
  const unsigned int q) const
{
  AssertIndexRange(q, this->n_quadrature_points);
  Assert(this->quadrature_points != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_quadrature_points"));
 
  // Cartesian/affine mesh: only first vertex of cell is stored, we must
  // compute it through the Jacobian (which is stored in non-inverted and
  // non-transposed form as index '1' in the jacobian field)
  if (is_face == false &&
      this->cell_type <= internal::MatrixFreeFunctions::affine)
    {
      Assert(this->jacobian != nullptr, ExcNotInitialized());
      Point<dim, Number> point = this->quadrature_points[0];
 
      const Tensor<2, dim, Number> &jac = this->jacobian[1];
      if (this->cell_type == internal::MatrixFreeFunctions::cartesian)
        for (unsigned int d = 0; d < dim; ++d)
          point[d] += jac[d][d] * static_cast<typename Number::value_type>(
                                    this->descriptor->quadrature.point(q)[d]);
      else
        for (unsigned int d = 0; d < dim; ++d)
          for (unsigned int e = 0; e < dim; ++e)
            point[d] += jac[d][e] * static_cast<typename Number::value_type>(
                                      this->descriptor->quadrature.point(q)[e]);
      return point;
    }
  else
    return this->quadrature_points[q];
}
 
 
 
template <int dim, typename Number, bool is_face>
inline DEAL_II_ALWAYS_INLINE Tensor<2, dim, Number>
FEEvaluationData<dim, Number, is_face>::inverse_jacobian(
  const unsigned int q_point) const
{
  AssertIndexRange(q_point, n_quadrature_points);
  Assert(jacobian != nullptr,
         internal::ExcMatrixFreeAccessToUninitializedMappingField(
           "update_gradients"));
  if (cell_type <= internal::MatrixFreeFunctions::affine)
    return jacobian[0];
  else
    return jacobian[q_point];
}
 
 
 
template <int dim, typename Number, bool is_face>
inline const Number *
FEEvaluationData<dim, Number, is_face>::begin_dof_values() const
{
  return values_dofs;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline Number *
FEEvaluationData<dim, Number, is_face>::begin_dof_values()
{
  if constexpr (running_in_debug_mode())
    {
      dof_values_initialized = true;
    }
  return values_dofs;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline const Number *
FEEvaluationData<dim, Number, is_face>::begin_values() const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(values_quad_initialized || values_quad_submitted,
             ExcNotInitialized());
    }
  return values_quad;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline Number *
FEEvaluationData<dim, Number, is_face>::begin_values()
{
  if constexpr (running_in_debug_mode())
    {
      values_quad_initialized = true;
      values_quad_submitted   = true;
    }
  return values_quad;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline const Number *
FEEvaluationData<dim, Number, is_face>::begin_gradients() const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(gradients_quad_initialized || gradients_quad_submitted,
             ExcNotInitialized());
    }
  return gradients_quad;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline Number *
FEEvaluationData<dim, Number, is_face>::begin_gradients()
{
  if constexpr (running_in_debug_mode())
    {
      gradients_quad_submitted   = true;
      gradients_quad_initialized = true;
    }
  return gradients_quad;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline const Number *
FEEvaluationData<dim, Number, is_face>::begin_hessians() const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(hessians_quad_initialized, ExcNotInitialized());
    }
  return hessians_quad;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline Number *
FEEvaluationData<dim, Number, is_face>::begin_hessians()
{
  if constexpr (running_in_debug_mode())
    {
      hessians_quad_initialized = true;
    }
  return hessians_quad;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline unsigned int
FEEvaluationData<dim, Number, is_face>::get_mapping_data_index_offset() const
{
  Assert(mapping_data != nullptr, ExcInternalError());
 
  if (dof_info == nullptr)
    return 0;
  else
    {
      AssertIndexRange(cell, mapping_data->data_index_offsets.size());
      return mapping_data->data_index_offsets[cell];
    }
}
 
 
 
template <int dim, typename Number, bool is_face>
inline internal::MatrixFreeFunctions::GeometryType
FEEvaluationData<dim, Number, is_face>::get_cell_type() const
{
  if constexpr (running_in_debug_mode())
    {
      Assert(is_reinitialized, ExcNotInitialized());
    }
  return cell_type;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline const typename FEEvaluationData<dim, Number, is_face>::ShapeInfoType &
FEEvaluationData<dim, Number, is_face>::get_shape_info() const
{
  Assert(data != nullptr, ExcInternalError());
  return *data;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline const internal::MatrixFreeFunctions::DoFInfo &
FEEvaluationData<dim, Number, is_face>::get_dof_info() const
{
  Assert(dof_info != nullptr,
         ExcMessage(
           "FEEvaluation was not initialized with a MatrixFree object!"));
  return *dof_info;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline const std::vector<unsigned int> &
FEEvaluationData<dim, Number, is_face>::get_internal_dof_numbering() const
{
  return data->lexicographic_numbering;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline unsigned int
FEEvaluationData<dim, Number, is_face>::get_quadrature_index() const
{
  return quadrature_index;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline unsigned int
FEEvaluationData<dim, Number, is_face>::get_current_cell_index() const
{
  if (is_face && dof_access_index ==
                   internal::MatrixFreeFunctions::DoFInfo::dof_access_cell)
    return cell * ReferenceCells::max_n_faces<dim>() + face_numbers[0];
  else
    return cell;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline unsigned int
FEEvaluationData<dim, Number, is_face>::get_active_fe_index() const
{
  return active_fe_index;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline unsigned int
FEEvaluationData<dim, Number, is_face>::get_active_quadrature_index() const
{
  return active_quad_index;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline unsigned int
FEEvaluationData<dim, Number, is_face>::get_first_selected_component() const
{
  return first_selected_component;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline ArrayView<Number>
FEEvaluationData<dim, Number, is_face>::get_scratch_data() const
{
  return scratch_data;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline std::uint8_t
FEEvaluationData<dim, Number, is_face>::get_face_no(const unsigned int v) const
{
  Assert(is_face, ExcNotInitialized());
  Assert(v == 0 || (cell != numbers::invalid_unsigned_int &&
                    dof_access_index ==
                      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
                    is_interior_face() == false),
         ExcMessage("All face numbers can only be queried for ECL at exterior "
                    "faces. Use get_face_no() in other cases."));
 
  return face_numbers[v];
}
 
 
 
template <int dim, typename Number, bool is_face>
inline unsigned int
FEEvaluationData<dim, Number, is_face>::get_subface_index() const
{
  return subface_index;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline std::uint8_t
FEEvaluationData<dim, Number, is_face>::get_face_orientation(
  const unsigned int v) const
{
  Assert(is_face, ExcNotInitialized());
  Assert(v == 0 || (cell != numbers::invalid_unsigned_int &&
                    dof_access_index ==
                      internal::MatrixFreeFunctions::DoFInfo::dof_access_cell &&
                    is_interior_face() == false),
         ExcMessage("All face numbers can only be queried for ECL at exterior "
                    "faces. Use get_face_no() in other cases."));
 
  return face_orientations[v];
}
 
 
 
template <int dim, typename Number, bool is_face>
inline internal::MatrixFreeFunctions::DoFInfo::DoFAccessIndex
FEEvaluationData<dim, Number, is_face>::get_dof_access_index() const
{
  return dof_access_index;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline bool
FEEvaluationData<dim, Number, is_face>::is_interior_face() const
{
  return interior_face;
}
 
 
 
template <int dim, typename Number, bool is_face>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FEEvaluationData<dim, Number, is_face>::quadrature_point_indices() const
{
  return std_cxx20::ranges::iota_view<unsigned int, unsigned int>(
    0U, n_quadrature_points);
}
 
 
 
namespace internal
{
  template <std::size_t N,
            typename VectorOfArrayType,
            typename ArrayType,
            typename FU>
  inline void
  process_cell_or_face_data(const std::array<unsigned int, N> indices,
                            VectorOfArrayType                &array,
                            ArrayType                        &out,
                            const FU                         &fu)
  {
    for (unsigned int i = 0; i < N; ++i)
      if (indices[i] != numbers::invalid_unsigned_int)
        {
          AssertIndexRange(indices[i] / N, array.size());
          fu(out[i], array[indices[i] / N][indices[i] % N]);
        }
  }
 
  template <std::size_t N, typename VectorOfArrayType, typename ArrayType>
  inline void
  set_valid_element_to_array(const std::array<unsigned int, N> indices,
                             const VectorOfArrayType          &array,
                             ArrayType                        &out)
  {
    AssertDimension(indices.size(), out.size());
    AssertDimension(indices.size(), array[0].size());
    // set all entries in array to a valid element
    std::size_t index = 0;
    for (; index < N; ++index)
      if (indices[index] != numbers::invalid_unsigned_int)
        break;
    for (unsigned int i = 0; i < N; ++i)
      out[i] = array[indices[index] / N][indices[index] % N];
  }
} // namespace internal
 
 
 
template <int dim, typename Number, bool is_face>
template <typename T>
inline T
FEEvaluationData<dim, Number, is_face>::read_cell_data(
  const AlignedVector<T> &array) const
{
  T out;
  internal::set_valid_element_to_array(this->get_cell_ids(), array, out);
  internal::process_cell_or_face_data(this->get_cell_ids(),
                                      array,
                                      out,
                                      [](auto &local, const auto &global) {
                                        local = global;
                                      });
  return out;
}
 
 
template <int dim, typename Number, bool is_face>
template <typename T, int D>
void
FEEvaluationData<dim, Number, is_face>::read_cell_data(
  const Table<D, T> &src,
  Table<D - 1, T>   &dst) const
{
  static_assert(D >= 2, "Table dimension must be at least 2");
 
 
  for (unsigned int d = 0; d < D - 1; ++d)
    Assert(dst.size(d) == src.size(d + 1),
           ExcMessage("Dimension mismatch between src and dst tables."));
 
  for (unsigned int lane = 0; lane < n_lanes; ++lane)
    if (this->get_cell_ids()[lane] != numbers::invalid_unsigned_int)
      {
        const unsigned int src_index  = this->get_cell_ids()[lane] / n_lanes;
        const unsigned int array_lane = this->get_cell_ids()[lane] % n_lanes;
 
        AssertIndexRange(src_index, src.size(0));
 
        // Copy data following the pattern:
        // dst[i][j][lane] = src[batch][i][j][lane]
        auto copy_recursively = [&](auto         self,
                                    auto       &&dst_slice,
                                    const auto  &src_slice,
                                    unsigned int depth) -> void {
          if constexpr (std::is_same_v<std::decay_t<decltype(dst_slice)>, T>)
            {
              // Base case: we've reached the innermost level (the array)
              dst_slice[lane] = src_slice[array_lane];
            }
          else
            {
              // Recursive case: iterate through this dimension. 'depth'
              // corresponds to the dimension index in the original 'src'
              // table.
              AssertIndexRange(depth, D);
              for (unsigned int i = 0; i < src.size(depth); ++i)
                self(self, dst_slice[i], src_slice[i], depth + 1);
            }
        };
 
        copy_recursively(copy_recursively, dst, src[src_index], 1);
      }
}
 
 
 
template <int dim, typename Number, bool is_face>
template <typename T>
inline void
FEEvaluationData<dim, Number, is_face>::set_cell_data(AlignedVector<T> &array,
                                                      const T &in) const
{
  internal::process_cell_or_face_data(this->get_cell_ids(),
                                      array,
                                      in,
                                      [](const auto &local, auto &global) {
                                        global = local;
                                      });
}
 
 
 
template <int dim, typename Number, bool is_face>
template <typename T>
inline T
FEEvaluationData<dim, Number, is_face>::read_face_data(
  const AlignedVector<T> &array) const
{
  T out;
  internal::set_valid_element_to_array(this->get_cell_ids(), array, out);
  internal::process_cell_or_face_data(this->get_face_ids(),
                                      array,
                                      out,
                                      [](auto &local, const auto &global) {
                                        local = global;
                                      });
  return out;
}
 
 
 
template <int dim, typename Number, bool is_face>
template <typename T>
inline void
FEEvaluationData<dim, Number, is_face>::set_face_data(AlignedVector<T> &array,
                                                      const T &in) const
{
  internal::process_cell_or_face_data(this->get_face_ids(),
                                      array,
                                      in,
                                      [](const auto &local, auto &global) {
                                        global = local;
                                      });
}
 
 
 
#endif // ifndef DOXYGEN
 
 
DEAL_II_NAMESPACE_CLOSE
 
#endif