field_container.h

Go to the documentation of this file.

// SPDX-FileCopyrightText: © 2025 PRISMS Center at the University of Michigan
// SPDX-License-Identifier: GNU Lesser General Public Version 2.1
 
#pragma once
 
#include <deal.II/base/config.h>
#include <deal.II/base/exceptions.h>
#include <deal.II/base/point.h>
#include <deal.II/base/symmetric_tensor.h>
#include <deal.II/base/tensor.h>
#include <deal.II/lac/vector.h>
#include <deal.II/matrix_free/evaluation_flags.h>
#include <deal.II/matrix_free/fe_evaluation.h>
#include <deal.II/matrix_free/matrix_free.h>
 
#include <prismspf/core/dependencies.h>
#include <prismspf/core/field_attributes.h>
#include <prismspf/core/solution_indexer.h>
#include <prismspf/core/solve_block.h>
#include <prismspf/core/system_wide.h>
#include <prismspf/core/type_enums.h>
#include <prismspf/core/types.h>
 
#include <prismspf/config.h>
 
#include <utility>
 
PRISMS_PF_BEGIN_NAMESPACE

template <unsigned int dim, unsigned int degree, typename number>
class FieldContainer
{
public:
  using ScalarValue = dealii::VectorizedArray<number>;

  using VectorValue = dealii::Tensor<1, dim, ScalarValue>;
 
  template <TensorRank Rank>
  using Value = std::conditional_t<Rank == TensorRank::Scalar || dim == 1,
                                   ScalarValue,
                                   dealii::Tensor<int(Rank), dim, ScalarValue>>;
 
  template <TensorRank Rank>
  using Gradient = dealii::Tensor<int(Rank) + 1, dim, ScalarValue>;
 
  template <TensorRank Rank>
  using Hessian = dealii::Tensor<int(Rank) + 2, dim, ScalarValue>;
 
  template <TensorRank Rank>
  using FEEval =
    dealii::FEEvaluation<dim,
                         degree,
                         degree + 1,
                         dealii::Tensor<int(Rank), dim>::n_independent_components,
                         number,
                         ScalarValue>;

  template <typename ValType>
  static constexpr TensorRank RankFromVal = []() constexpr
  {
    if constexpr (std::is_same_v<ValType, ScalarValue> ||
                  std::is_same_v<ValType, typename ScalarValue::value_type>)
      {
        return TensorRank::Scalar;
      }
    else
      {
        return TensorRank(ValType::rank);
      }
  }();

  template <typename GradType>
  static constexpr TensorRank RankFromGrad = []() constexpr
  {
    if constexpr (std::is_same_v<GradType, ScalarValue>)
      {
        return TensorRank::Scalar;
      }
    else
      {
        return TensorRank(GradType::rank - 1);
      }
  }();

  template <TensorRank Rank>
  struct FEEValuationDeps
  {
    using FEEDepPair = std::pair<FEEval<Rank>, EvalFlags>;

    using FEEDepPairPtr = std::shared_ptr<FEEDepPair>;

    FEEDepPairPtr fe_eval;

    FEEDepPairPtr fe_eval_src_dst;

    std::vector<FEEDepPairPtr> fe_eval_old;

    EvalFlags integration_flags = EvalFlags::nothing;

    const SolutionLevel<dim, number> *solution_level = nullptr;

    unsigned int field_index = -1;

    unsigned int block_index = -1;

    FEEValuationDeps() = default;
 
    FEEValuationDeps(const MatrixFree<dim, number>    &matrix_free,
                     const Dependency                 &dependency,
                     const SolutionLevel<dim, number> &_solution_level,
                     unsigned int                      _field_index,
                     unsigned int                      _block_index,
                     bool                              is_dst);
 
    template <DependencyType type>
    const FEEval<Rank> &
    get() const;
 
    template <DependencyType type>
    FEEval<Rank> &
    get();
 
    const FEEval<Rank> &
    get(DependencyType type) const;
 
    FEEval<Rank> &
    get(DependencyType type);
 
    void
    reinit(unsigned int cell);
 
    void
    eval(const BlockVector<number> *_src_solutions, bool plain);
 
    void
    reinit_and_eval(unsigned int               cell,
                    const BlockVector<number> *_src_solutions,
                    bool                       plain);
 
    void
    eval_without_read();
 
    void
    integrate();
 
    void
    distribute(BlockVector<number> *dst_solutions);
 
    void
    integrate_and_distribute(BlockVector<number> *dst_solutions);
  };

  FieldContainer(const std::vector<FieldAttributes> &_field_attributes,
                 const SolutionIndexer<dim, number> &_solution_indexer,
                 unsigned int                        _relative_level,
                 const DependencyMap                &dependency_map,
                 const SolveBlock                   &_solve_block,
                 const MatrixFree<dim, number>      &matrix_free);

  void
  reinit(unsigned int cell);

  void
  eval(const BlockVector<number> *src_solutions, bool plain);

  void
  reinit_and_eval(unsigned int               cell,
                  const BlockVector<number> *src_solutions,
                  bool                       plain);

  void
  eval_without_read();

  void
  integrate();

  void
  distribute(BlockVector<number> *dst_solutions);

  void
  distribute(unsigned int field_index, BlockVector<number> *dst_solutions);

  void
  integrate_and_distribute(BlockVector<number> *dst_solutions);

  void
  set_q_point(unsigned int q);

  template <TensorRank Rank, DependencyType type>
  [[nodiscard]] Value<Rank>
  get_value(Types::Index field_index) const;

  template <TensorRank Rank>
  [[nodiscard]] Value<Rank>
  get_value(Types::Index field_index, DependencyType type) const;

  template <TensorRank Rank, DependencyType type>
  [[nodiscard]] Gradient<Rank>
  get_gradient(Types::Index field_index) const;

  template <TensorRank Rank>
  [[nodiscard]] Gradient<Rank>
  get_gradient(Types::Index field_index, DependencyType type) const;

  template <TensorRank Rank, DependencyType type>
  [[nodiscard]] Hessian<Rank>
  get_hessian(Types::Index field_index) const;

  template <TensorRank Rank>
  [[nodiscard]] Hessian<Rank>
  get_hessian(Types::Index field_index, DependencyType type) const;

  template <TensorRank Rank, DependencyType type>
  [[nodiscard]] Gradient<Rank>
  get_hessian_diagonal(Types::Index field_index) const;

  template <TensorRank Rank>
  [[nodiscard]] Gradient<Rank>
  get_hessian_diagonal(Types::Index field_index, DependencyType type) const;

  template <TensorRank Rank, DependencyType type>
  [[nodiscard]] Value<Rank>
  get_laplacian(Types::Index field_index) const;

  template <TensorRank Rank>
  [[nodiscard]] Value<Rank>
  get_laplacian(Types::Index field_index, DependencyType type) const;

  template <TensorRank Rank, DependencyType type>
  [[nodiscard]] ScalarValue
  get_divergence(Types::Index field_index) const;

  template <TensorRank Rank>
  [[nodiscard]] ScalarValue
  get_divergence(Types::Index field_index, DependencyType type) const;

  template <TensorRank Rank, DependencyType type>
  [[nodiscard]] dealii::SymmetricTensor<2, dim, ScalarValue>
  get_symmetric_gradient(Types::Index field_index) const;

  template <TensorRank Rank>
  [[nodiscard]] dealii::SymmetricTensor<2, dim, ScalarValue>
  get_symmetric_gradient(Types::Index field_index, DependencyType type) const;

  template <TensorRank Rank, DependencyType type>
  [[nodiscard]] dealii::Tensor<1, (dim == 2 ? 1 : dim), ScalarValue>
  get_curl(Types::Index field_index) const;

  template <TensorRank Rank>
  [[nodiscard]] dealii::Tensor<1, (dim == 2 ? 1 : dim), ScalarValue>
  get_curl(Types::Index field_index, DependencyType type) const;

  [[nodiscard]] dealii::Point<dim, ScalarValue>
  get_q_point_location() const;

  [[nodiscard]] ScalarValue
  get_element_volume() const;

  [[nodiscard]] unsigned int
  get_n_q_points() const;

  template <typename ValType>
  void
  set_value_term(Types::Index field_index, const ValType &val);

  template <typename GradType>
  void
  set_gradient_term(Types::Index field_index, const GradType &val);

  template <typename ValType>
  void
  get_dof_value_to(ValType     &destination,
                   Types::Index field_index,
                   unsigned int dof_index);

  template <typename ValType>
  void
  submit_dof_value(Types::Index field_index, const ValType &val, unsigned int dof_index);
 
  constexpr static unsigned int dofs_per_component =
    FEEval<Scalar>::static_dofs_per_component;
 
private:
  template <TensorRank Rank>
  std::vector<FEEValuationDeps<Rank>> &
  get_relevant_feeval_vector();
 
  template <TensorRank Rank>
  const std::vector<FEEValuationDeps<Rank>> &
  get_relevant_feeval_vector() const;

  void
  feevaluation_size_valid(Types::Index field_index) const;

  void
  feevaluation_exists(Types::Index field_index, Types::Index dependency_index) const;

  void
  access_valid(Types::Index                             field_index,
               DependencyType                           dependency_type,
               dealii::EvaluationFlags::EvaluationFlags flag) const;

  void
  submission_valid(Types::Index field_index, DependencyType dependency_type) const;

  const std::vector<FieldAttributes> *field_attributes_ptr;

  const SolutionIndexer<dim, number> *solution_indexer;

  std::vector<FEEValuationDeps<TensorRank::Scalar>> feeval_deps_scalar;

  std::vector<FEEValuationDeps<TensorRank::Vector>> feeval_deps_vector;

  const SolveBlock *solve_block;

  FEEval<TensorRank::Scalar> shared_feeval_scalar;

  unsigned int relative_level;

  unsigned int q_point = 0;
};
 
inline static const std::map<DependencyType, std::string> dependency_type_to_string = {
  {DependencyType::Current,  "current"},
  {DependencyType::OldOne,   "old_1"  },
  {DependencyType::OldTwo,   "old_2"  },
  {DependencyType::OldThree, "old_3"  },
  {DependencyType::OldFour,  "old_4"  },
  {DependencyType::SRC,      "src"    },
  {DependencyType::DST,      "dst"    }
};
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::FEEValuationDeps(
  const MatrixFree<dim, number>    &matrix_free,
  const Dependency                 &dependency,
  const SolutionLevel<dim, number> &_solution_level,
  unsigned int                      _field_index,
  unsigned int                      _block_index,
  bool                              is_dst)
  : solution_level(&_solution_level)
  , field_index(_field_index)
  , block_index(_block_index)
{
  // Make an FEEvaluation if the current solution needs to be evaluated
  if (dependency.flag)
    {
      fe_eval = std::make_shared<FEEDepPair>(FEEval<Rank>(matrix_free, field_index),
                                             dependency.flag);
    }
  // Make FEEvaluations for the the old solutions
  fe_eval_old.resize(dependency.old_flags.size(), nullptr);
  for (unsigned int age = 0; age < dependency.old_flags.size(); ++age)
    {
      if (dependency.old_flags.at(age)) // kinda redundant... maybe remove?
        {
          fe_eval_old[age] =
            std::make_shared<FEEDepPair>(FEEval<Rank>(matrix_free, field_index),
                                         dependency.old_flags.at(age));
        }
    }
  if (dependency.src_flag || is_dst)
    {
      fe_eval_src_dst =
        std::make_shared<FEEDepPair>(FEEval<Rank>(matrix_free, field_index),
                                     dependency.src_flag);
    }
}
 
class ExcDepNotInitialized : public dealii::ExceptionBase
{
public:
  explicit ExcDepNotInitialized(DependencyType _dependency_type)
    : dependency_type(_dependency_type)
  {}
 
  const char *
  what() const noexcept override
  {
    return "Access was attempted for a field that was not declared as a "
           "dependency for the current solve block.\n";
  }
 
  void
  print_info(std::ostream &out) const override
  {
    out << what() << std::flush;
  }
 
private:
  DependencyType dependency_type;
};
#ifdef DEBUG
#  define AssertAccessible(fe_eval_pair_ptr, dependency_type)                   \
    AssertThrowDebug(block_index != -1, ExcDepNotInitialized(dependency_type)); \
    AssertThrowDebug(fe_eval_pair_ptr, ExcDepNotInitialized(dependency_type));
#else
#  define AssertAccessible(fe_eval_pair_ptr, dependency_type)
#endif
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
template <DependencyType Type>
inline DEAL_II_ALWAYS_INLINE const typename FieldContainer<dim, degree, number>::
  template FEEval<Rank> &
  FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::get() const
{
  if constexpr (Type == DependencyType::SRC || Type == DependencyType::DST)
    {
      AssertAccessible(fe_eval_src_dst, Type);
      return fe_eval_src_dst->first;
    }
  else if constexpr (Type == DependencyType::Current)
    {
      AssertAccessible(fe_eval, Type);
      return fe_eval->first;
    }
  else
    {
      AssertThrowDebug(int(Type) - 1 < fe_eval_old.size(), ExcDepNotInitialized(Type));
      AssertAccessible(fe_eval_old[int(Type) - 1], Type);
      return fe_eval_old[int(Type) - 1]->first;
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
template <DependencyType Type>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template FEEval<Rank> &
  FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::get()
{
  if constexpr (Type == DependencyType::SRC || Type == DependencyType::DST)
    {
      AssertAccessible(fe_eval_src_dst, Type);
      return fe_eval_src_dst->first;
    }
  else if constexpr (Type == DependencyType::Current)
    {
      AssertAccessible(fe_eval, Type);
      return fe_eval->first;
    }
  else
    {
      AssertThrowDebug(int(Type) - 1 < fe_eval_old.size(), ExcDepNotInitialized(Type));
      AssertAccessible(fe_eval_old[int(Type) - 1], Type);
      return fe_eval_old[int(Type) - 1]->first;
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE const typename FieldContainer<dim, degree, number>::
  template FEEval<Rank> &
  FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::get(
    DependencyType type) const
{
  if (type == DependencyType::SRC || type == DependencyType::DST)
    {
      AssertAccessible(fe_eval_src_dst, type);
      return fe_eval_src_dst->first;
    }
  if (type == DependencyType::Current)
    {
      AssertAccessible(fe_eval, type);
      return fe_eval->first;
    }
  {
    AssertThrowDebug(int(type) - 1 < fe_eval_old.size(), ExcDepNotInitialized(type));
    AssertAccessible(fe_eval_old[int(type) - 1], type);
    return fe_eval_old[type - 1]->first;
  }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template FEEval<Rank> &
  FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::get(DependencyType type)
{
  if (type == DependencyType::SRC || type == DependencyType::DST)
    {
      AssertAccessible(fe_eval_src_dst, type);
      return fe_eval_src_dst->first;
    }
  if (type == DependencyType::Current)
    {
      AssertAccessible(fe_eval, type);
      return fe_eval->first;
    }
  {
    AssertThrowDebug(int(type) - 1 < fe_eval_old.size(), ExcDepNotInitialized(type));
    AssertAccessible(fe_eval_old[int(type) - 1], type);
    return fe_eval_old[int(type) - 1]->first;
  }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline void
FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::reinit(unsigned int cell)
{
  if (fe_eval)
    {
      fe_eval->first.reinit(cell);
    }
  for (auto &old_fe_eval : fe_eval_old)
    {
      if (old_fe_eval)
        {
          old_fe_eval->first.reinit(cell);
        }
    }
  if (fe_eval_src_dst)
    {
      fe_eval_src_dst->first.reinit(cell);
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline void
FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::eval(
  const BlockVector<number> *_src_solutions,
  bool                       plain)
{
  // NOTE: `read_dof_values_plain` must be called here so that constraints aren't
  // implicitly applied. This allows us to have inhomogeneous constraints.
  if (fe_eval)
    {
      fe_eval->first.read_dof_values_plain(solution_level->solutions.block(block_index));
      fe_eval->first.evaluate(fe_eval->second);
    }
  for (unsigned int age = 0; age < fe_eval_old.size(); ++age)
    {
      if (FEEDepPairPtr &old_fe_eval = fe_eval_old[age])
        {
          old_fe_eval->first.read_dof_values_plain(
            solution_level->old_solutions[age].block(block_index));
          old_fe_eval->first.evaluate(old_fe_eval->second);
        }
    }
  if (fe_eval_src_dst && fe_eval_src_dst->second != EvalFlags::nothing)
    {
      if (plain)
        {
          fe_eval_src_dst->first.read_dof_values_plain(
            _src_solutions->block(block_index));
        }
      else
        {
          fe_eval_src_dst->first.read_dof_values(_src_solutions->block(block_index));
        }
      fe_eval_src_dst->first.evaluate(fe_eval_src_dst->second);
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline void
FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::reinit_and_eval(
  unsigned int               cell,
  const BlockVector<number> *_src_solutions,
  bool                       plain)
{
  // NOTE: `read_dof_values_plain` must be called here so that constraints aren't
  // implicitly applied. This allows us to have inhomogeneous constraints.
  if (fe_eval)
    {
      fe_eval->first.reinit(cell);
      fe_eval->first.read_dof_values_plain(solution_level->solutions.block(block_index));
      fe_eval->first.evaluate(fe_eval->second);
    }
  for (unsigned int age = 0; age < fe_eval_old.size(); ++age)
    {
      if (FEEDepPairPtr &old_fe_eval = fe_eval_old[age])
        {
          old_fe_eval->first.reinit(cell);
          old_fe_eval->first.read_dof_values_plain(
            solution_level->old_solutions[age].block(block_index));
          old_fe_eval->first.evaluate(old_fe_eval->second);
        }
    }
  if (fe_eval_src_dst)
    {
      fe_eval_src_dst->first.reinit(cell);
      if (fe_eval_src_dst->second != EvalFlags::nothing)
        {
          if (plain)
            {
              fe_eval_src_dst->first.read_dof_values_plain(
                _src_solutions->block(block_index));
            }
          else
            {
              fe_eval_src_dst->first.read_dof_values(_src_solutions->block(block_index));
            }
          fe_eval_src_dst->first.evaluate(fe_eval_src_dst->second);
        }
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline void
FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::eval_without_read()
{
  if (fe_eval_src_dst && fe_eval_src_dst->second != EvalFlags::nothing)
    {
      fe_eval_src_dst->first.evaluate(fe_eval_src_dst->second);
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline void
FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::integrate()
{
  if (fe_eval_src_dst)
    {
      fe_eval_src_dst->first.integrate(integration_flags);
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline void
FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::distribute(
  BlockVector<number> *dst_solutions)
{
  if (fe_eval_src_dst)
    {
      fe_eval_src_dst->first.distribute_local_to_global(
        dst_solutions->block(block_index));
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline void
FieldContainer<dim, degree, number>::FEEValuationDeps<Rank>::integrate_and_distribute(
  BlockVector<number> *dst_solutions)
{
  if (fe_eval_src_dst)
    {
      fe_eval_src_dst->first.integrate_scatter(integration_flags,
                                               dst_solutions->block(block_index));
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline void
FieldContainer<dim, degree, number>::reinit(unsigned int cell)
{
  for (auto &fe_eval : feeval_deps_scalar)
    {
      fe_eval.reinit(cell);
    }
  for (auto &fe_eval : feeval_deps_vector)
    {
      fe_eval.reinit(cell);
    }
  shared_feeval_scalar.reinit(cell);
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline void
FieldContainer<dim, degree, number>::eval(const BlockVector<number> *src_solutions,
                                          bool                       plain)
{
  for (auto &fe_eval : feeval_deps_scalar)
    {
      fe_eval.eval(src_solutions, plain);
    }
  for (auto &fe_eval : feeval_deps_vector)
    {
      fe_eval.eval(src_solutions, plain);
    }
  // Don't eval `shared_feeval_scalar` because we only use it for information.
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline void
FieldContainer<dim, degree, number>::reinit_and_eval(
  unsigned int               cell,
  const BlockVector<number> *src_solutions,
  bool                       plain)
{
  for (auto &fe_eval : feeval_deps_scalar)
    {
      fe_eval.reinit_and_eval(cell, src_solutions, plain);
    }
  for (auto &fe_eval : feeval_deps_vector)
    {
      fe_eval.reinit_and_eval(cell, src_solutions, plain);
    }
  // Don't eval `shared_feeval_scalar` because we only use it for information.
  shared_feeval_scalar.reinit(cell);
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline void
FieldContainer<dim, degree, number>::eval_without_read()
{
  for (auto &fe_eval : feeval_deps_scalar)
    {
      fe_eval.eval_without_read();
    }
  for (auto &fe_eval : feeval_deps_vector)
    {
      fe_eval.eval_without_read();
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline void
FieldContainer<dim, degree, number>::integrate()
{
  const std::vector<FieldAttributes> &field_attributes = *field_attributes_ptr;
  for (const Types::Index &field_index : solve_block->field_indices)
    {
      if (field_attributes[field_index].field_type == TensorRank::Scalar)
        {
          feeval_deps_scalar[field_index].integrate();
        }
      else /* vector */
        {
          feeval_deps_vector[field_index].integrate();
        }
    }
  // Don't integrate `shared_feeval_scalar` because we only use it for information.
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline void
FieldContainer<dim, degree, number>::distribute(BlockVector<number> *dst_solutions)
{
  const std::vector<FieldAttributes> &field_attributes = *field_attributes_ptr;
  for (const Types::Index &field_index : solve_block->field_indices)
    {
      if (field_attributes[field_index].field_type == TensorRank::Scalar)
        {
          feeval_deps_scalar[field_index].distribute(dst_solutions);
        }
      else /* vector */
        {
          feeval_deps_vector[field_index].distribute(dst_solutions);
        }
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline void
FieldContainer<dim, degree, number>::distribute(unsigned int         field_index,
                                                BlockVector<number> *dst_solutions)
{
  const std::vector<FieldAttributes> &field_attributes = *field_attributes_ptr;
  if (field_attributes[field_index].field_type == TensorRank::Scalar)
    {
      feeval_deps_scalar[field_index].distribute(dst_solutions);
    }
  else /* vector */
    {
      feeval_deps_vector[field_index].distribute(dst_solutions);
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline void
FieldContainer<dim, degree, number>::integrate_and_distribute(
  BlockVector<number> *dst_solutions)
{
  const std::vector<FieldAttributes> &field_attributes = *field_attributes_ptr;
  for (const Types::Index &field_index : solve_block->field_indices)
    {
      if (field_attributes[field_index].field_type == TensorRank::Scalar)
        {
          feeval_deps_scalar[field_index].integrate_and_distribute(dst_solutions);
        }
      else /* vector */
        {
          feeval_deps_vector[field_index].integrate_and_distribute(dst_solutions);
        }
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline DEAL_II_ALWAYS_INLINE void
FieldContainer<dim, degree, number>::set_q_point(unsigned int q)
{
  q_point = q;
}
 
// there are two catches we can do here.
// 1. Dependencies for the dependency type (current, old, src/dst) don't exist.
// 2. Dependency is not initialized for values/gradients.
// We catch these separately to give more informative error messages.
#define GetterTempl(dependency_type) template get<dependency_type>()
#define GetterNoTempl(dependency_type) get(dependency_type)
#define ReturnGetter(get_handle, Rank, field_index, dependency_type, getter)        \
  try                                                                               \
    {                                                                               \
      AssertThrowDebug((field_index) < get_relevant_feeval_vector<Rank>().size(),   \
                       dealii::ExcMessage("Error: Field index " +                   \
                                          std::to_string(field_index) +             \
                                          " is not associated with any field."));   \
      return get_relevant_feeval_vector<Rank>()[field_index]                        \
        .getter(dependency_type)                                                    \
        .get_handle(q_point);                                                       \
    }                                                                               \
  catch (const ExcDepNotInitialized &e)                                             \
    {                                                                               \
      std::cerr << "Error when trying to access field with index " << (field_index) \
                << " and dependency type "                                          \
                << dependency_type_to_string.at(dependency_type) << ":\n"           \
                << e.what()                                                         \
                << "Ensure that each dependency is requested in the solve block.\n" \
                << std::flush;                                                      \
      throw;                                                                        \
    }                                                                               \
  catch (const dealii::internal::ExcAccessToUninitializedField &e)                  \
    {                                                                               \
      std::cerr << "Error when trying to access field with index " << (field_index) \
                << " and dependency type "                                          \
                << dependency_type_to_string.at(dependency_type)                    \
                << ":\nAccess was attempted for values or gradients that were "     \
                   "not requested.\n"                                               \
                << std::flush;                                                      \
      throw;                                                                        \
    }
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank, DependencyType type>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template Value<Rank>
  FieldContainer<dim, degree, number>::get_value(Types::Index field_index) const
{
  ReturnGetter(get_value, Rank, field_index, type, GetterTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template Value<Rank>
  FieldContainer<dim, degree, number>::get_value(Types::Index   field_index,
                                                 DependencyType type) const
{
  ReturnGetter(get_value, Rank, field_index, type, GetterNoTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank, DependencyType type>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template Gradient<Rank>
  FieldContainer<dim, degree, number>::get_gradient(Types::Index field_index) const
{
  ReturnGetter(get_gradient, Rank, field_index, type, GetterTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template Gradient<Rank>
  FieldContainer<dim, degree, number>::get_gradient(Types::Index   field_index,
                                                    DependencyType type) const
{
  ReturnGetter(get_gradient, Rank, field_index, type, GetterNoTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank, DependencyType type>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template Hessian<Rank>
  FieldContainer<dim, degree, number>::get_hessian(Types::Index field_index) const
{
  ReturnGetter(get_hessian, Rank, field_index, type, GetterTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template Hessian<Rank>
  FieldContainer<dim, degree, number>::get_hessian(Types::Index   field_index,
                                                   DependencyType type) const
{
  ReturnGetter(get_hessian, Rank, field_index, type, GetterNoTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank, DependencyType type>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template Gradient<Rank>
  FieldContainer<dim, degree, number>::get_hessian_diagonal(
    Types::Index field_index) const
{
  ReturnGetter(get_hessian_diagonal, Rank, field_index, type, GetterTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template Gradient<Rank>
  FieldContainer<dim, degree, number>::get_hessian_diagonal(Types::Index   field_index,
                                                            DependencyType type) const
{
  ReturnGetter(get_hessian_diagonal, Rank, field_index, type, GetterNoTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank, DependencyType type>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template Value<Rank>
  FieldContainer<dim, degree, number>::get_laplacian(Types::Index field_index) const
{
  ReturnGetter(get_laplacian, Rank, field_index, type, GetterTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE
  typename FieldContainer<dim, degree, number>::template Value<Rank>
  FieldContainer<dim, degree, number>::get_laplacian(Types::Index   field_index,
                                                     DependencyType type) const
{
  ReturnGetter(get_laplacian, Rank, field_index, type, GetterNoTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank, DependencyType type>
inline DEAL_II_ALWAYS_INLINE typename FieldContainer<dim, degree, number>::ScalarValue
FieldContainer<dim, degree, number>::get_divergence(Types::Index field_index) const
{
  static_assert(Rank == 1, "Divergences are only available for vector fields");
  ReturnGetter(get_divergence, Rank, field_index, type, GetterTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE typename FieldContainer<dim, degree, number>::ScalarValue
FieldContainer<dim, degree, number>::get_divergence(Types::Index   field_index,
                                                    DependencyType type) const
{
  static_assert(Rank == 1, "Divergences are only available for vector fields");
  ReturnGetter(get_divergence, Rank, field_index, type, GetterNoTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank, DependencyType type>
inline DEAL_II_ALWAYS_INLINE dealii::
  SymmetricTensor<2, dim, typename FieldContainer<dim, degree, number>::ScalarValue>
  FieldContainer<dim, degree, number>::get_symmetric_gradient(
    Types::Index field_index) const
{
  static_assert(Rank == 1, "Symmetric gradients are only available for vector fields");
  ReturnGetter(get_symmetric_gradient, Rank, field_index, type, GetterTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE dealii::
  SymmetricTensor<2, dim, typename FieldContainer<dim, degree, number>::ScalarValue>
  FieldContainer<dim, degree, number>::get_symmetric_gradient(Types::Index   field_index,
                                                              DependencyType type) const
{
  static_assert(Rank == 1, "Symmetric gradients are only available for vector fields");
  ReturnGetter(get_symmetric_gradient, Rank, field_index, type, GetterNoTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank, DependencyType type>
inline DEAL_II_ALWAYS_INLINE
  dealii::Tensor<1,
                 (dim == 2 ? 1 : dim),
                 typename FieldContainer<dim, degree, number>::ScalarValue>
  FieldContainer<dim, degree, number>::get_curl(Types::Index field_index) const
{
  static_assert(Rank == 1, "Curl is only available for vector fields");
  AssertThrowDebug(
    dim > 1,
    dealii::ExcMessage(
      "Curl is only available for vector fields with dimension greater than 1."));
  ReturnGetter(get_curl, Rank, field_index, type, GetterTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE
  dealii::Tensor<1,
                 (dim == 2 ? 1 : dim),
                 typename FieldContainer<dim, degree, number>::ScalarValue>
  FieldContainer<dim, degree, number>::get_curl(Types::Index   field_index,
                                                DependencyType type) const
{
  static_assert(Rank == 1, "Curl is only available for vector fields");
  AssertThrowDebug(
    dim > 1,
    dealii::ExcMessage(
      "Curl is only available for vector fields with dimension greater than 1."));
  ReturnGetter(get_curl, Rank, field_index, type, GetterNoTempl);
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline DEAL_II_ALWAYS_INLINE
  dealii::Point<dim, typename FieldContainer<dim, degree, number>::ScalarValue>
  FieldContainer<dim, degree, number>::get_q_point_location() const
{
  return shared_feeval_scalar.quadrature_point(q_point);
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline DEAL_II_ALWAYS_INLINE typename FieldContainer<dim, degree, number>::ScalarValue
FieldContainer<dim, degree, number>::get_element_volume() const
{
  return shared_feeval_scalar.JxW(q_point) /
         SystemWide<dim, degree>::quadrature.weight(q_point);
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline DEAL_II_ALWAYS_INLINE unsigned int
FieldContainer<dim, degree, number>::get_n_q_points() const
{
  return shared_feeval_scalar.n_q_points;
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <typename ValType>
inline DEAL_II_ALWAYS_INLINE void
FieldContainer<dim, degree, number>::set_value_term(Types::Index   field_index,
                                                    const ValType &val)
{
  AssertThrowDebug((field_index) <
                     get_relevant_feeval_vector<RankFromVal<ValType>>().size(),
                   dealii::ExcMessage("Error: Field index " +
                                      std::to_string(field_index) +
                                      " is not associated with any field."));
  auto &relevant_feeval_vector =
    get_relevant_feeval_vector<RankFromVal<ValType>>()[field_index];
  try
    {
      relevant_feeval_vector.template get<DependencyType::DST>().submit_value(val,
                                                                              q_point);
      relevant_feeval_vector.integration_flags |= dealii::EvaluationFlags::values;
    }
  catch (...)
    {
      std::cerr << "Error when trying to submit value for field with index "
                << (field_index)
                << ": Error: Submission for field not part of this solve block.\n"
                << std::flush;
      throw;
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <typename GradType>
inline DEAL_II_ALWAYS_INLINE void
FieldContainer<dim, degree, number>::set_gradient_term(Types::Index    field_index,
                                                       const GradType &val)
{
  AssertThrowDebug((field_index) <
                     get_relevant_feeval_vector<RankFromGrad<GradType>>().size(),
                   dealii::ExcMessage("Error: Field index " +
                                      std::to_string(field_index) +
                                      " is not associated with any field."));
  auto &relevant_feeval_vector =
    get_relevant_feeval_vector<RankFromGrad<GradType>>()[field_index];
  try
    {
      relevant_feeval_vector.template get<DependencyType::DST>().submit_gradient(val,
                                                                                 q_point);
      relevant_feeval_vector.integration_flags |= dealii::EvaluationFlags::gradients;
    }
  catch (...)
    {
      std::cerr << "Error when trying to submit gradient for field with index "
                << (field_index)
                << ": Error: Submission for field not part of this solve block.\n"
                << std::flush;
      throw;
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <typename ValType>
inline DEAL_II_ALWAYS_INLINE void
FieldContainer<dim, degree, number>::get_dof_value_to(ValType     &destination,
                                                      Types::Index field_index,
                                                      unsigned int dof_index)
{
  AssertThrowDebug((field_index) <
                     get_relevant_feeval_vector<RankFromVal<ValType>>().size(),
                   dealii::ExcMessage("Error: Field index " +
                                      std::to_string(field_index) +
                                      " is not associated with any field."));
  auto &relevant_feeval_vector =
    get_relevant_feeval_vector<RankFromVal<ValType>>()[field_index];
  try
    {
      destination =
        relevant_feeval_vector.template get<DependencyType::DST>().get_dof_value(
          dof_index);
    }
  catch (...)
    {
      std::cerr << "Error when trying to get dof value for field with index "
                << (field_index) << ": Error: Field not part of this solve block.\n"
                << std::flush;
      throw;
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <typename ValType>
inline DEAL_II_ALWAYS_INLINE void
FieldContainer<dim, degree, number>::submit_dof_value(Types::Index   field_index,
                                                      const ValType &val,
                                                      unsigned int   dof_index)
{
  AssertThrowDebug((field_index) <
                     get_relevant_feeval_vector<RankFromVal<ValType>>().size(),
                   dealii::ExcMessage("Error: Field index " +
                                      std::to_string(field_index) +
                                      " is not associated with any field."));
  auto &relevant_feevals =
    get_relevant_feeval_vector<RankFromVal<ValType>>()[field_index];
  try
    {
      relevant_feevals.template get<DependencyType::DST>().submit_dof_value(val,
                                                                            dof_index);
    }
  catch (...)
    {
      std::cerr << "Error when trying to submit dof value for field with index "
                << (field_index)
                << ": Error: Submission for field not part of this solve block.\n"
                << std::flush;
      throw;
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE std::vector<
  typename FieldContainer<dim, degree, number>::template FEEValuationDeps<Rank>> &
FieldContainer<dim, degree, number>::get_relevant_feeval_vector()
{
  if constexpr (Rank == TensorRank::Scalar || dim == 1)
    {
      return feeval_deps_scalar;
    }
  else if constexpr (Rank == TensorRank::Vector)
    {
      return feeval_deps_vector;
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
template <TensorRank Rank>
inline DEAL_II_ALWAYS_INLINE const std::vector<
  typename FieldContainer<dim, degree, number>::template FEEValuationDeps<Rank>> &
FieldContainer<dim, degree, number>::get_relevant_feeval_vector() const
{
  if constexpr (Rank == TensorRank::Scalar || dim == 1)
    {
      return feeval_deps_scalar;
    }
  else if constexpr (Rank == TensorRank::Vector)
    {
      return feeval_deps_vector;
    }
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline DEAL_II_ALWAYS_INLINE void
FieldContainer<dim, degree, number>::feevaluation_size_valid(
  [[maybe_unused]] Types::Index field_index) const
{
  // TODO
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline DEAL_II_ALWAYS_INLINE void
FieldContainer<dim, degree, number>::feevaluation_exists(
  [[maybe_unused]] Types::Index field_index,
  [[maybe_unused]] Types::Index dependency_index) const
{
  // TODO
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline DEAL_II_ALWAYS_INLINE void
FieldContainer<dim, degree, number>::access_valid(
  [[maybe_unused]] Types::Index                             field_index,
  [[maybe_unused]] DependencyType                           dependency_type,
  [[maybe_unused]] dealii::EvaluationFlags::EvaluationFlags flag) const
{
  // TODO
}
 
template <unsigned int dim, unsigned int degree, typename number>
inline DEAL_II_ALWAYS_INLINE void
FieldContainer<dim, degree, number>::submission_valid(
  [[maybe_unused]] Types::Index   field_index,
  [[maybe_unused]] DependencyType dependency_type) const
{
  // TODO
}
 
#undef AssertAccessible
#undef ReturnGetter
#undef GetterTempl
#undef GetterNoTempl
 
PRISMS_PF_END_NAMESPACE