ad_helpers.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2016 - 2017 by the deal.II authors
//
// This file is part of the deal.II library.
//
// The deal.II library is free software; you can use it, redistribute
// it, and/or modify it under the terms of the GNU Lesser General
// Public License as published by the Free Software Foundation; either
// version 2.1 of the License, or (at your option) any later version.
// The full text of the license can be found in the file LICENSE at
// the top level of the deal.II distribution.
//
// ---------------------------------------------------------------------

#ifndef dealii_differentiation_ad_ad_helpers_h
#define dealii_differentiation_ad_ad_helpers_h

#include <deal.II/base/config.h>

#if defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)

#  include <deal.II/base/numbers.h>
#  include <deal.II/base/symmetric_tensor.h>
#  include <deal.II/base/tensor.h>

#  include <deal.II/differentiation/ad/ad_drivers.h>
#  include <deal.II/differentiation/ad/ad_number_traits.h>
#  include <deal.II/differentiation/ad/adolc_number_types.h>
#  include <deal.II/differentiation/ad/adolc_product_types.h>
#  include <deal.II/differentiation/ad/sacado_number_types.h>
#  include <deal.II/differentiation/ad/sacado_product_types.h>

#  include <deal.II/fe/fe_values_extractors.h>

#  include <deal.II/lac/full_matrix.h>
#  include <deal.II/lac/vector.h>

#  include <algorithm>
#  include <iostream>
#  include <iterator>
#  include <numeric>
#  include <set>

DEAL_II_NAMESPACE_OPEN

namespace Differentiation
{
  namespace AD
  {
    template <enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType = double>
    class HelperBase
    {
    public:
      using scalar_type =
        typename AD::NumberTraits<ScalarType, ADNumberTypeCode>::scalar_type;

      using ad_type =
        typename AD::NumberTraits<ScalarType, ADNumberTypeCode>::ad_type;


      HelperBase(const unsigned int n_independent_variables,
                 const unsigned int n_dependent_variables);

      virtual ~HelperBase() = default;



      std::size_t
      n_independent_variables() const;

      std::size_t
      n_dependent_variables() const;

      void
      print(std::ostream &stream) const;

      void
      print_values(std::ostream &stream) const;

      void
      print_tape_stats(const typename Types<ad_type>::tape_index tape_index,
                       std::ostream &                            stream) const;



      static void
      configure_tapeless_mode(const unsigned int n_independent_variables,
                              const bool ensure_persistent_setting = true);



      virtual void
      reset(const unsigned int n_independent_variables =
              ::numbers::invalid_unsigned_int,
            const unsigned int n_dependent_variables =
              ::numbers::invalid_unsigned_int,
            const bool clear_registered_tapes = true);

      bool
      is_recording() const;

      typename Types<ad_type>::tape_index
      active_tape_index() const;

      bool
      is_registered_tape(
        const typename Types<ad_type>::tape_index tape_index) const;

      void
      set_tape_buffer_sizes(
        const typename Types<ad_type>::tape_buffer_sizes obufsize = 64 * 1024 *
                                                                    1024,
        const typename Types<ad_type>::tape_buffer_sizes lbufsize = 64 * 1024 *
                                                                    1024,
        const typename Types<ad_type>::tape_buffer_sizes vbufsize = 64 * 1024 *
                                                                    1024,
        const typename Types<ad_type>::tape_buffer_sizes tbufsize = 64 * 1024 *
                                                                    1024);

      bool
      start_recording_operations(
        const typename Types<ad_type>::tape_index tape_index,
        const bool                                overwrite_tape = false,
        const bool keep_independent_values                       = true);

      void
      stop_recording_operations(const bool write_tapes_to_file = false);

      void
      activate_recorded_tape(
        const typename Types<ad_type>::tape_index tape_index);

      bool
      recorded_tape_requires_retaping(
        const typename Types<ad_type>::tape_index tape_index) const;

      bool
      active_tape_requires_retaping() const;

      void
      clear_active_tape();


    protected:

      TapedDrivers<ad_type, scalar_type> taped_driver;

      TapelessDrivers<ad_type, scalar_type> tapeless_driver;

      void
      activate_tape(const typename Types<ad_type>::tape_index tape_index,
                    const bool                                read_mode);



      mutable std::vector<scalar_type> independent_variable_values;

      mutable std::vector<ad_type> independent_variables;

      std::vector<bool> registered_independent_variable_values;

      mutable std::vector<bool> registered_marked_independent_variables;

      void
      reset_registered_independent_variables();

      void
      set_sensitivity_value(const unsigned int index, const scalar_type &value);

      void
      mark_independent_variable(const unsigned int index, ad_type &out) const;

      void
      finalize_sensitive_independent_variables() const;

      void
      initialize_non_sensitive_independent_variable(const unsigned int index,
                                                    ad_type &out) const;

      unsigned int
      n_registered_independent_variables() const;



      std::vector<ad_type> dependent_variables;

      std::vector<bool> registered_marked_dependent_variables;

      void
      reset_registered_dependent_variables(const bool flag = false);

      unsigned int
      n_registered_dependent_variables() const;

      void
      register_dependent_variable(const unsigned int index,
                                  const ad_type &    func);


    }; // class HelperBase



    template <enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType = double>
    class CellLevelBase : public HelperBase<ADNumberTypeCode, ScalarType>
    {
    public:
      using scalar_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::scalar_type;

      using ad_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::ad_type;


      CellLevelBase(const unsigned int n_independent_variables,
                    const unsigned int n_dependent_variables);

      virtual ~CellLevelBase() = default;



      void
      register_dof_values(const std::vector<scalar_type> &dof_values);

      template <typename VectorType>
      void
      register_dof_values(
        const VectorType &                                  values,
        const std::vector<::types::global_dof_index> &local_dof_indices);

      const std::vector<ad_type> &
      get_sensitive_dof_values() const;



      void
      set_dof_values(const std::vector<scalar_type> &dof_values);

      template <typename VectorType>
      void
      set_dof_values(
        const VectorType &                                  values,
        const std::vector<::types::global_dof_index> &local_dof_indices);



      virtual void
      compute_residual(Vector<scalar_type> &residual) const = 0;

      virtual void
      compute_linearization(FullMatrix<scalar_type> &linearization) const = 0;


    }; // class CellLevelBase



    template <enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType = double>
    class EnergyFunctional : public CellLevelBase<ADNumberTypeCode, ScalarType>
    {
    public:
      using scalar_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::scalar_type;

      using ad_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::ad_type;


      EnergyFunctional(const unsigned int n_independent_variables);

      virtual ~EnergyFunctional() = default;



      void
      register_energy_functional(const ad_type &energy);

      scalar_type
      compute_energy() const;

      void
      compute_residual(Vector<scalar_type> &residual) const override;

      virtual void
      compute_linearization(
        FullMatrix<scalar_type> &linearization) const override;


    }; // class EnergyFunctional



    template <enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType = double>
    class ResidualLinearization
      : public CellLevelBase<ADNumberTypeCode, ScalarType>
    {
    public:
      using scalar_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::scalar_type;

      using ad_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::ad_type;


      ResidualLinearization(const unsigned int n_independent_variables,
                            const unsigned int n_dependent_variables);

      virtual ~ResidualLinearization() = default;



      void
      register_residual_vector(const std::vector<ad_type> &residual);

      virtual void
      compute_residual(Vector<scalar_type> &residual) const override;

      virtual void
      compute_linearization(
        FullMatrix<scalar_type> &linearization) const override;


    }; // class ResidualLinearization



    namespace internal
    {
      template <int dim, typename ExtractorType>
      struct Extractor;


      template <int dim>
      struct Extractor<dim, FEValuesExtractors::Scalar>
      {
        static const unsigned int n_components = 1;

        static const unsigned int rank = 0;

        template <typename NumberType>
        using tensor_type = Tensor<rank, dim, NumberType>;

        static_assert(
          n_components == tensor_type<double>::n_independent_components,
          "The number of components doesn't match that of the corresponding tensor type.");
        static_assert(
          rank == tensor_type<double>::rank,
          "The rank doesn't match that of the corresponding tensor type.");

        // Note: FEValuesViews::Scalar::tensor_type is double, so we can't
        // use it (FEValuesViews) in this context.
        // In fact, sadly, all FEValuesViews objects expect doubles as value
        // types.
        template <typename NumberType>
        using value_type = NumberType;

        template <typename NumberType>
        using gradient_type = Tensor<rank + 1, dim, NumberType>; // NumberType;

        static inline unsigned int
        first_component(const FEValuesExtractors::Scalar &extractor)
        {
          return extractor.component;
        }

        static bool
        symmetric_component(const unsigned int unrolled_index)
        {
          (void)unrolled_index;
          return false;
        }

        template <typename IndexType = unsigned int, int rank_in>
        static IndexType
        local_component(const TableIndices<rank_in> &table_indices,
                        const unsigned int           column_offset)
        {
          Assert(column_offset <= rank_in, ExcInternalError());
          (void)table_indices;
          (void)column_offset;
          return 0;
        }
      };


      template <int dim>
      struct Extractor<dim, FEValuesExtractors::Vector>
      {
        static const unsigned int n_components = dim;

        static const unsigned int rank = 1;

        template <typename NumberType>
        using tensor_type = Tensor<rank, dim, NumberType>;

        static_assert(
          n_components == tensor_type<double>::n_independent_components,
          "The number of components doesn't match that of the corresponding tensor type.");
        static_assert(
          rank == tensor_type<double>::rank,
          "The rank doesn't match that of the corresponding tensor type.");

        template <typename NumberType>
        using value_type = tensor_type<NumberType>;

        template <typename NumberType>
        using gradient_type = Tensor<rank + 1, dim, NumberType>;

        static inline unsigned int
        first_component(const FEValuesExtractors::Vector &extractor)
        {
          return extractor.first_vector_component;
        }

        static bool
        symmetric_component(const unsigned int unrolled_index)
        {
          (void)unrolled_index;
          return false;
        }

        template <int rank_in>
        static TableIndices<rank>
        table_index_view(const TableIndices<rank_in> &table_indices,
                         const unsigned int           column_offset)
        {
          Assert(0 + column_offset < rank_in, ExcInternalError());
          return TableIndices<rank>(table_indices[column_offset]);
        }

        template <typename IndexType = unsigned int, int rank_in>
        static IndexType
        local_component(const TableIndices<rank_in> &table_indices,
                        const unsigned int           column_offset)
        {
          static_assert(
            rank_in >= rank,
            "Cannot extract more table indices than the input table has!");
          using TensorType = tensor_type<double>;
          return TensorType::component_to_unrolled_index(
            table_index_view(table_indices, column_offset));
        }
      };


      template <int dim>
      struct Extractor<dim, FEValuesExtractors::Tensor<1>>
      {
        static const unsigned int n_components =
          Tensor<1, dim>::n_independent_components;

        static const unsigned int rank = 1;

        template <typename NumberType>
        using tensor_type = Tensor<rank, dim, NumberType>;

        template <typename NumberType>
        using value_type = tensor_type<NumberType>;

        template <typename NumberType>
        using gradient_type = Tensor<rank + 1, dim, NumberType>;

        static inline unsigned int
        first_component(const FEValuesExtractors::Tensor<1> &extractor)
        {
          return extractor.first_tensor_component;
        }

        static bool
        symmetric_component(const unsigned int unrolled_index)
        {
          (void)unrolled_index;
          return false;
        }

        template <int rank_in>
        static TableIndices<rank>
        table_index_view(const TableIndices<rank_in> &table_indices,
                         const unsigned int           column_offset)
        {
          Assert(column_offset < rank_in, ExcInternalError());
          return TableIndices<rank>(table_indices[column_offset]);
        }

        template <typename IndexType = unsigned int, int rank_in>
        static IndexType
        local_component(const TableIndices<rank_in> &table_indices,
                        const unsigned int           column_offset)
        {
          static_assert(
            rank_in >= rank,
            "Cannot extract more table indices than the input table has!");
          using TensorType = tensor_type<double>;
          return TensorType::component_to_unrolled_index(
            table_index_view(table_indices, column_offset));
        }
      };


      template <int dim>
      struct Extractor<dim, FEValuesExtractors::Tensor<2>>
      {
        static const unsigned int n_components =
          Tensor<2, dim>::n_independent_components;

        static const unsigned int rank = Tensor<2, dim>::rank;

        template <typename NumberType>
        using tensor_type = Tensor<rank, dim, NumberType>;

        template <typename NumberType>
        using value_type = tensor_type<NumberType>;

        template <typename NumberType>
        using gradient_type = Tensor<rank + 1, dim, NumberType>;

        static inline unsigned int
        first_component(const FEValuesExtractors::Tensor<2> &extractor)
        {
          return extractor.first_tensor_component;
        }

        static bool
        symmetric_component(const unsigned int unrolled_index)
        {
          (void)unrolled_index;
          return false;
        }

        template <int rank_in>
        static TableIndices<rank>
        table_index_view(const TableIndices<rank_in> &table_indices,
                         const unsigned int           column_offset)
        {
          Assert(column_offset < rank_in, ExcInternalError());
          Assert(column_offset + 1 < rank_in, ExcInternalError());
          return TableIndices<rank>(table_indices[column_offset],
                                    table_indices[column_offset + 1]);
        }

        template <typename IndexType = unsigned int, int rank_in>
        static IndexType
        local_component(const TableIndices<rank_in> &table_indices,
                        const unsigned int           column_offset)
        {
          static_assert(
            rank_in >= rank,
            "Cannot extract more table indices than the input table has!");
          using TensorType = tensor_type<double>;
          return TensorType::component_to_unrolled_index(
            table_index_view(table_indices, column_offset));
        }
      };


      template <int dim>
      struct Extractor<dim, FEValuesExtractors::SymmetricTensor<2>>
      {
        static const unsigned int n_components =
          SymmetricTensor<2, dim>::n_independent_components;

        static const unsigned int rank = SymmetricTensor<2, dim>::rank;

        template <typename NumberType>
        using tensor_type = SymmetricTensor<rank, dim, NumberType>;

        template <typename NumberType>
        using value_type = tensor_type<NumberType>;

        template <typename NumberType>
        using gradient_type = Tensor<rank + 1, dim, NumberType>;

        static inline unsigned int
        first_component(const FEValuesExtractors::SymmetricTensor<2> &extractor)
        {
          return extractor.first_tensor_component;
        }

        static bool
        symmetric_component(const unsigned int unrolled_index)
        {
          const TableIndices<2> table_indices =
            tensor_type<double>::unrolled_to_component_indices(unrolled_index);
          return table_indices[0] != table_indices[1];
        }

        template <int rank_in>
        static TableIndices<rank>
        table_index_view(const TableIndices<rank_in> &table_indices,
                         const unsigned int           column_offset)
        {
          Assert(column_offset < rank_in, ExcInternalError());
          Assert(column_offset + 1 < rank_in, ExcInternalError());
          return TableIndices<rank>(table_indices[column_offset],
                                    table_indices[column_offset + 1]);
        }

        template <typename IndexType = unsigned int, int rank_in>
        static IndexType
        local_component(const TableIndices<rank_in> &table_indices,
                        const unsigned int           column_offset)
        {
          static_assert(
            rank_in >= rank,
            "Cannot extract more table indices than the input table has!");
          using TensorType = tensor_type<double>;
          return TensorType::component_to_unrolled_index(
            table_index_view(table_indices, column_offset));
        }
      };


      template <int dim, typename NumberType, typename ExtractorType>
      struct ScalarFieldGradient
      {
        using type =
          typename Extractor<dim,
                             ExtractorType>::template tensor_type<NumberType>;
      };


      template <typename ExtractorType_Row, typename ExtractorType_Col>
      struct HessianType
      {
        template <int rank, int dim, typename NumberType>
        using type = Tensor<rank, dim, NumberType>;
      };


      template <>
      struct HessianType<FEValuesExtractors::SymmetricTensor<2>,
                         FEValuesExtractors::Scalar>
      {
        template <int /*rank*/, int dim, typename NumberType>
        using type = SymmetricTensor<2 /*rank*/, dim, NumberType>;
      };


      template <>
      struct HessianType<FEValuesExtractors::Scalar,
                         FEValuesExtractors::SymmetricTensor<2>>
      {
        template <int /*rank*/, int dim, typename NumberType>
        using type = SymmetricTensor<2 /*rank*/, dim, NumberType>;
      };


      template <>
      struct HessianType<FEValuesExtractors::SymmetricTensor<2>,
                         FEValuesExtractors::SymmetricTensor<2>>
      {
        template <int /*rank*/, int dim, typename NumberType>
        using type = SymmetricTensor<4 /*rank*/, dim, NumberType>;
      };


      template <int dim,
                typename NumberType,
                typename ExtractorType_Row,
                typename ExtractorType_Col>
      struct ScalarFieldHessian
      {
        static const int rank = Extractor<dim, ExtractorType_Row>::rank +
                                Extractor<dim, ExtractorType_Col>::rank;

        using type =
          typename HessianType<ExtractorType_Row, ExtractorType_Col>::
            template type<rank, dim, NumberType>;
      };


      template <int dim, typename NumberType, typename ExtractorType_Field>
      using VectorFieldValue =
        ScalarFieldGradient<dim, NumberType, ExtractorType_Field>;


      template <int dim,
                typename NumberType,
                typename ExtractorType_Field,
                typename ExtractorType_Derivative>
      using VectorFieldJacobian = ScalarFieldHessian<dim,
                                                     NumberType,
                                                     ExtractorType_Field,
                                                     ExtractorType_Derivative>;


      template <int dim,
                typename IndexType = unsigned int,
                typename ExtractorType>
      std::vector<IndexType>
      extract_field_component_indices(const ExtractorType &extractor,
                                      const bool ignore_symmetries = true)
      {
        (void)ignore_symmetries;
        const IndexType n_components =
          internal::Extractor<dim, ExtractorType>::n_components;
        const IndexType comp_first =
          internal::Extractor<dim, ExtractorType>::first_component(extractor);
        std::vector<IndexType> indices(n_components);
        std::iota(indices.begin(), indices.end(), comp_first);
        return indices;
      }


      template <int dim, typename IndexType = unsigned int>
      std::vector<IndexType>
      extract_field_component_indices(
        const FEValuesExtractors::SymmetricTensor<2> &extractor_symm_tensor,
        const bool                                    ignore_symmetries = true)
      {
        using ExtractorType = FEValuesExtractors::SymmetricTensor<2>;
        const IndexType n_components =
          internal::Extractor<dim, ExtractorType>::n_components;
        if (ignore_symmetries == true)
          {
            const IndexType comp_first =
              internal::Extractor<dim, ExtractorType>::first_component(
                extractor_symm_tensor);
            std::vector<IndexType> indices(n_components);
            std::iota(indices.begin(), indices.end(), comp_first);
            return indices;
          }
        else
          {
            // First get all of the indices of the non-symmetric tensor
            const FEValuesExtractors::Tensor<2> extractor_tensor(
              extractor_symm_tensor.first_tensor_component);
            std::vector<IndexType> indices =
              extract_field_component_indices<dim>(extractor_tensor, true);

            // Then we overwrite any illegal entries with the equivalent indices
            // from the symmetric tensor
            for (unsigned int i = 0; i < indices.size(); ++i)
              {
                if (indices[i] >= n_components)
                  {
                    const TableIndices<2> ti_tensor =
                      Tensor<2, dim>::unrolled_to_component_indices(indices[i]);
                    const IndexType sti_new_index =
                      SymmetricTensor<2, dim>::component_to_unrolled_index(
                        ti_tensor);
                    indices[i] = sti_new_index;
                  }
              }

            return indices;
          }
      }


      template <typename TensorType, typename NumberType>
      inline void
      set_tensor_entry(TensorType &        t,
                       const unsigned int &unrolled_index,
                       const NumberType &  value)
      {
        // Where possible, set values using TableIndices
        Assert(unrolled_index < t.n_independent_components,
               ExcIndexRange(unrolled_index, 0, t.n_independent_components));
        t[TensorType::unrolled_to_component_indices(unrolled_index)] = value;
      }


      template <int dim, typename NumberType>
      inline void set_tensor_entry(Tensor<0, dim, NumberType> &t,
                                   const unsigned int &        unrolled_index,
                                   const NumberType &          value)
      {
        Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1));
        (void)unrolled_index;
        t = value;
      }


      template <typename NumberType>
      inline void
      set_tensor_entry(NumberType &        t,
                       const unsigned int &unrolled_index,
                       const NumberType &  value)
      {
        Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1));
        (void)unrolled_index;
        t = value;
      }


      template <int dim, typename NumberType>
      inline void set_tensor_entry(SymmetricTensor<4, dim, NumberType> &t,
                                   const unsigned int &unrolled_index_row,
                                   const unsigned int &unrolled_index_col,
                                   const NumberType &  value)
      {
        // Fourth order symmetric tensors require a specialized interface
        // to extract values.
        using SubTensorType = SymmetricTensor<2, dim, NumberType>;
        Assert(unrolled_index_row < SubTensorType::n_independent_components,
               ExcIndexRange(unrolled_index_row,
                             0,
                             SubTensorType::n_independent_components));
        Assert(unrolled_index_col < SubTensorType::n_independent_components,
               ExcIndexRange(unrolled_index_col,
                             0,
                             SubTensorType::n_independent_components));
        const TableIndices<2> indices_row =
          SubTensorType::unrolled_to_component_indices(unrolled_index_row);
        const TableIndices<2> indices_col =
          SubTensorType::unrolled_to_component_indices(unrolled_index_col);
        t[indices_row[0]][indices_row[1]][indices_col[0]][indices_col[1]] =
          value;
      }


      template <int rank,
                int dim,
                typename NumberType,
                template <int, int, typename> class TensorType>
      inline NumberType
      get_tensor_entry(const TensorType<rank, dim, NumberType> &t,
                       const unsigned int &                     unrolled_index)
      {
        // Where possible, get values using TableIndices
        Assert(unrolled_index < t.n_independent_components,
               ExcIndexRange(unrolled_index, 0, t.n_independent_components));
        return t[TensorType<rank, dim, NumberType>::
                   unrolled_to_component_indices(unrolled_index)];
      }


      template <int dim,
                typename NumberType,
                template <int, int, typename> class TensorType>
      inline NumberType
      get_tensor_entry(const TensorType<0, dim, NumberType> &t,
                       const unsigned int &                  unrolled_index)
      {
        Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1));
        (void)unrolled_index;
        return t;
      }


      template <typename NumberType>
      inline const NumberType &
      get_tensor_entry(const NumberType &t, const unsigned int &unrolled_index)
      {
        Assert(unrolled_index == 0, ExcIndexRange(unrolled_index, 0, 1));
        (void)unrolled_index;
        return t;
      }


      template <int rank,
                int dim,
                typename NumberType,
                template <int, int, typename> class TensorType>
      inline NumberType &
      get_tensor_entry(TensorType<rank, dim, NumberType> &t,
                       const unsigned int &               unrolled_index)
      {
        // Where possible, get values using TableIndices
        Assert(unrolled_index < t.n_independent_components,
               ExcIndexRange(unrolled_index, 0, t.n_independent_components));
        return t[TensorType<rank, dim, NumberType>::
                   unrolled_to_component_indices(unrolled_index)];
      }


      template <int dim,
                typename NumberType,
                template <int, int, typename> class TensorType>
      NumberType &get_tensor_entry(TensorType<0, dim, NumberType> &t,
                                   const unsigned int &            index)
      {
        Assert(index == 0, ExcIndexRange(index, 0, 1));
        (void)index;
        return t;
      }


      template <typename NumberType>
      inline NumberType &
      get_tensor_entry(NumberType &t, const unsigned int &index)
      {
        Assert(index == 0, ExcIndexRange(index, 0, 1));
        (void)index;
        return t;
      }

    } // namespace internal



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType = double>
    class PointLevelFunctionsBase
      : public HelperBase<ADNumberTypeCode, ScalarType>
    {
    public:
      static const unsigned int dimension = dim;

      using scalar_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::scalar_type;

      using ad_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::ad_type;


      PointLevelFunctionsBase(const unsigned int n_independent_variables,
                              const unsigned int n_dependent_variables);

      virtual ~PointLevelFunctionsBase() = default;



      virtual void
      reset(const unsigned int n_independent_variables =
              ::numbers::invalid_unsigned_int,
            const unsigned int n_dependent_variables =
              ::numbers::invalid_unsigned_int,
            const bool clear_registered_tapes = true) override;

      void
      register_independent_variables(const std::vector<scalar_type> &values);

      template <typename ValueType, typename ExtractorType>
      void
      register_independent_variable(const ValueType &    value,
                                    const ExtractorType &extractor);

      const std::vector<ad_type> &
      get_sensitive_variables() const;

      /*
       * Extract a subset of the independent variables as represented by
       * auto-differentiable numbers. These are the independent
       * variables @f$\mathbf{A} \subset \mathbf{X}@f$ at which the dependent values
       * are evaluated and differentiated.
       *
       * This function indicates to the AD library that implements the
       * auto-differentiable number type that operations performed on these
       * numbers are to be tracked so they are considered "sensitive"
       * variables. This is, therefore, the set of variables with which one
       * would then perform computations, and based on which one can then
       * extract both the value of the function and its derivatives with the
       * member functions below. The values of the components of the returned
       * object are initialized to the values set with
       * register_independent_variable().
       *
       * @param[in] extractor An extractor associated with the input field
       * variables. This effectively defines which components of the global set
       * of independent variables this field is associated with.
       * @return An object of auto-differentiable type numbers. The return type is
       * based on the input extractor, and will be either a scalar,
       * Tensor<1,dim>, Tensor<2,dim>, or SymmetricTensor<2,dim>.
       *
       * @note For taped AD numbers, this operation is only valid in recording mode.
       */
      template <typename ExtractorType>
      typename internal::Extractor<dim,
                                   ExtractorType>::template tensor_type<ad_type>
      get_sensitive_variables(const ExtractorType &extractor) const;



      void
      set_independent_variables(const std::vector<scalar_type> &values);

      template <typename ValueType, typename ExtractorType>
      void
      set_independent_variable(const ValueType &    value,
                               const ExtractorType &extractor);


    protected:

      void
      set_sensitivity_value(const unsigned int index,
                            const bool         symmetric_component,
                            const scalar_type &value);

      bool
      is_symmetric_independent_variable(const unsigned int index) const;

      unsigned int
      n_symmetric_independent_variables() const;


    private:

      std::vector<bool> symmetric_independent_variables;


    }; // class PointLevelFunctionsBase



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType = double>
    class ScalarFunction
      : public PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>
    {
    public:
      using scalar_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::scalar_type;

      using ad_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::ad_type;


      ScalarFunction(const unsigned int n_independent_variables);

      virtual ~ScalarFunction() = default;



      void
      register_dependent_variable(const ad_type &func);

      scalar_type
      compute_value() const;

      void
      compute_gradient(Vector<scalar_type> &gradient) const;

      void
      compute_hessian(FullMatrix<scalar_type> &hessian) const;

      template <typename ExtractorType_Row>
      static typename internal::
        ScalarFieldGradient<dim, scalar_type, ExtractorType_Row>::type
        extract_gradient_component(const Vector<scalar_type> &gradient,
                                   const ExtractorType_Row &  extractor_row);

      template <typename ExtractorType_Row, typename ExtractorType_Col>
      static typename internal::ScalarFieldHessian<dim,
                                                   scalar_type,
                                                   ExtractorType_Row,
                                                   ExtractorType_Col>::type
      extract_hessian_component(const FullMatrix<scalar_type> &hessian,
                                const ExtractorType_Row &      extractor_row,
                                const ExtractorType_Col &      extractor_col);

      static Tensor<0, dim, scalar_type>
      extract_hessian_component(
        const FullMatrix<scalar_type> &   hessian,
        const FEValuesExtractors::Scalar &extractor_row,
        const FEValuesExtractors::Scalar &extractor_col);

      static SymmetricTensor<4, dim, scalar_type>
      extract_hessian_component(
        const FullMatrix<scalar_type> &               hessian,
        const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
        const FEValuesExtractors::SymmetricTensor<2> &extractor_col);


    }; // class ScalarFunction



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType = double>
    class VectorFunction
      : public PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>
    {
    public:
      using scalar_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::scalar_type;

      using ad_type =
        typename HelperBase<ADNumberTypeCode, ScalarType>::ad_type;


      VectorFunction(const unsigned int n_independent_variables,
                     const unsigned int n_dependent_variables);

      virtual ~VectorFunction() = default;



      void
      register_dependent_variables(const std::vector<ad_type> &funcs);

      template <typename ValueType, typename ExtractorType>
      void
      register_dependent_variable(const ValueType &    funcs,
                                  const ExtractorType &extractor);

      void
      compute_values(Vector<scalar_type> &values) const;

      void
      compute_jacobian(FullMatrix<scalar_type> &jacobian) const;


      template <typename ExtractorType_Row>
      static typename internal::
        VectorFieldValue<dim, scalar_type, ExtractorType_Row>::type
        extract_value_component(const Vector<scalar_type> &values,
                                const ExtractorType_Row &  extractor_row);

      template <typename ExtractorType_Row, typename ExtractorType_Col>
      static typename internal::VectorFieldJacobian<dim,
                                                    scalar_type,
                                                    ExtractorType_Row,
                                                    ExtractorType_Col>::type
      extract_jacobian_component(const FullMatrix<scalar_type> &jacobian,
                                 const ExtractorType_Row &      extractor_row,
                                 const ExtractorType_Col &      extractor_col);

      static Tensor<0, dim, scalar_type>
      extract_jacobian_component(
        const FullMatrix<scalar_type> &   jacobian,
        const FEValuesExtractors::Scalar &extractor_row,
        const FEValuesExtractors::Scalar &extractor_col);

      static SymmetricTensor<4, dim, scalar_type>
      extract_jacobian_component(
        const FullMatrix<scalar_type> &               jacobian,
        const FEValuesExtractors::SymmetricTensor<2> &extractor_row,
        const FEValuesExtractors::SymmetricTensor<2> &extractor_col);


    }; // class VectorFunction


  } // namespace AD
} // namespace Differentiation


/* ----------------- inline and template functions ----------------- */


#  ifndef DOXYGEN

namespace Differentiation
{
  namespace AD
  {
    /* ----------------- CellLevelBase ----------------- */



    template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
    template <typename VectorType>
    void
    CellLevelBase<ADNumberTypeCode, ScalarType>::register_dof_values(
      const VectorType &                                  values,
      const std::vector<::types::global_dof_index> &local_dof_indices)
    {
      // This is actually the same thing the set_dof_values() function,
      // in the sense that we simply populate our array of independent values
      // with a meaningful number. However, in this case we need to double check
      // that we're not registering these variables twice
      Assert(
        local_dof_indices.size() == this->n_independent_variables(),
        ExcMessage(
          "Degree of freedom index vector size does not match number of independent variables"));
      for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
        {
          Assert(this->registered_independent_variable_values[i] == false,
                 ExcMessage("Independent variables already registered."));
        }
      set_dof_values(values, local_dof_indices);
    }



    template <enum AD::NumberTypes ADNumberTypeCode, typename ScalarType>
    template <typename VectorType>
    void
    CellLevelBase<ADNumberTypeCode, ScalarType>::set_dof_values(
      const VectorType &                                  values,
      const std::vector<::types::global_dof_index> &local_dof_indices)
    {
      Assert(local_dof_indices.size() == this->n_independent_variables(),
             ExcMessage(
               "Vector size does not match number of independent variables"));
      for (unsigned int i = 0; i < this->n_independent_variables(); ++i)
        HelperBase<ADNumberTypeCode, ScalarType>::set_sensitivity_value(
          i, values[local_dof_indices[i]]);
    }



    /* ----------------- PointLevelFunctionsBase ----------------- */



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType>
    template <typename ValueType, typename ExtractorType>
    void
    PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>::
      register_independent_variable(const ValueType &    value,
                                    const ExtractorType &extractor)
    {
      // This is actually the same thing as the set_independent_variable
      // function, in the sense that we simply populate our array of independent
      // values with a meaningful number. However, in this case we need to
      // double check that we're not registering these variables twice
#    ifdef DEBUG
      const std::vector<unsigned int> index_set(
        internal::extract_field_component_indices<dim>(extractor));
      for (unsigned int i = 0; i < index_set.size(); ++i)
        {
          Assert(
            this->registered_independent_variable_values[index_set[i]] == false,
            ExcMessage(
              "Overlapping indices for independent variables. "
              "One or more indices associated with the field that "
              "is being registered as an independent variable have "
              "already been associated with another field. This suggests "
              "that the component offsets used to construct their counterpart "
              "extractors are incompatible with one another. Make sure that "
              "the first component for each extractor properly takes into "
              "account the dimensionality of the preceding fields."));
        }
#    endif
      set_independent_variable(value, extractor);
    }



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType>
    template <typename ValueType, typename ExtractorType>
    void
    PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>::
      set_independent_variable(const ValueType &    value,
                               const ExtractorType &extractor)
    {
      const std::vector<unsigned int> index_set(
        internal::extract_field_component_indices<dim>(extractor));
      for (unsigned int i = 0; i < index_set.size(); ++i)
        {
          set_sensitivity_value(
            index_set[i],
            internal::Extractor<dim, ExtractorType>::symmetric_component(i),
            internal::get_tensor_entry(value, i));
        }
    }



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType>
    template <typename ExtractorType>
    typename internal::Extractor<dim, ExtractorType>::template tensor_type<
      typename HelperBase<ADNumberTypeCode, ScalarType>::ad_type>
    PointLevelFunctionsBase<dim, ADNumberTypeCode, ScalarType>::
      get_sensitive_variables(const ExtractorType &extractor) const
    {
      if (ADNumberTraits<ad_type>::is_taped == true)
        {
          Assert(this->active_tape_index() !=
                   Numbers<ad_type>::invalid_tape_index,
                 ExcMessage("Invalid tape index"));
        }

      // If necessary, finalize the internally stored vector of
      // AD numbers that represents the independent variables
      this->finalize_sensitive_independent_variables();
      Assert(this->independent_variables.size() ==
               this->n_independent_variables(),
             ExcDimensionMismatch(this->independent_variables.size(),
                                  this->n_independent_variables()));

      const std::vector<unsigned int> index_set(
        internal::extract_field_component_indices<dim>(extractor));
      typename internal::Extractor<dim,
                                   ExtractorType>::template tensor_type<ad_type>
        out;

      for (unsigned int i = 0; i < index_set.size(); ++i)
        {
          const unsigned int index = index_set[i];
          Assert(index < this->n_independent_variables(), ExcInternalError());
          Assert(this->registered_independent_variable_values[index] == true,
                 ExcInternalError());
          internal::get_tensor_entry(out, i) =
            this->independent_variables[index];
        }

      return out;
    }



    /* ----------------- ScalarFunction ----------------- */



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType>
    template <typename ExtractorType_Row>
    typename internal::ScalarFieldGradient<
      dim,
      typename ScalarFunction<dim, ADNumberTypeCode, ScalarType>::scalar_type,
      ExtractorType_Row>::type
    ScalarFunction<dim, ADNumberTypeCode, ScalarType>::
      extract_gradient_component(const Vector<scalar_type> &gradient,
                                 const ExtractorType_Row &  extractor_row)
    {
      // NOTE: The order of components must be consistently defined throughout
      // this class.
      typename internal::
        ScalarFieldGradient<dim, scalar_type, ExtractorType_Row>::type out;

      // Get indexsets for the subblock from which we wish to extract the
      // gradient values
      const std::vector<unsigned int> row_index_set(
        internal::extract_field_component_indices<dim>(extractor_row));
      Assert(out.n_independent_components == row_index_set.size(),
             ExcMessage("Not all tensor components have been extracted!"));
      for (unsigned int r = 0; r < row_index_set.size(); ++r)
        internal::set_tensor_entry(out, r, gradient[row_index_set[r]]);

      return out;
    }



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType>
    template <typename ExtractorType_Row, typename ExtractorType_Col>
    typename internal::ScalarFieldHessian<
      dim,
      typename ScalarFunction<dim, ADNumberTypeCode, ScalarType>::scalar_type,
      ExtractorType_Row,
      ExtractorType_Col>::type
    ScalarFunction<dim, ADNumberTypeCode, ScalarType>::
      extract_hessian_component(const FullMatrix<scalar_type> &hessian,
                                const ExtractorType_Row &      extractor_row,
                                const ExtractorType_Col &      extractor_col)
    {
      using InternalHessian      = internal::ScalarFieldHessian<dim,
                                                           scalar_type,
                                                           ExtractorType_Row,
                                                           ExtractorType_Col>;
      using InternalExtractorRow = internal::Extractor<dim, ExtractorType_Row>;
      using InternalExtractorCol = internal::Extractor<dim, ExtractorType_Col>;
      using HessianType          = typename InternalHessian::type;

      // NOTE: The order of components must be consistently defined throughout
      // this class.
      HessianType out;

      // Get indexsets for the subblocks from which we wish to extract the
      // Hessian values
      // NOTE: Here we have to do some clever accounting when the
      // one extractor is a symmetric Tensor and the other is not, e.g.
      // <SymmTensor,Vector>. In this scenario the return type is a
      // non-symmetric Tensor<3,dim> but we have to fetch information from a
      // SymmTensor row/column that has too few entries to fill the output
      // tensor. So we must duplicate the relevant entries in the row/column
      // indexset to fetch off-diagonal components that are Otherwise
      // non-existent in a SymmTensor.
      const std::vector<unsigned int> row_index_set(
        internal::extract_field_component_indices<dim>(
          extractor_row, false /*ignore_symmetries*/));
      const std::vector<unsigned int> col_index_set(
        internal::extract_field_component_indices<dim>(
          extractor_col, false /*ignore_symmetries*/));

      for (unsigned int index = 0;
           index < HessianType::n_independent_components;
           ++index)
        {
          const TableIndices<HessianType::rank> ti_out =
            HessianType::unrolled_to_component_indices(index);
          const unsigned int r =
            InternalExtractorRow::local_component(ti_out, 0);
          const unsigned int c =
            InternalExtractorCol::local_component(ti_out,
                                                  InternalExtractorRow::rank);

          internal::set_tensor_entry(
            out, index, hessian[row_index_set[r]][col_index_set[c]]);
        }

      return out;
    }



    /* ----------------- VectorFunction ----------------- */



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType>
    template <typename ValueType, typename ExtractorType>
    void
    VectorFunction<dim, ADNumberTypeCode, ScalarType>::
      register_dependent_variable(const ValueType &    funcs,
                                  const ExtractorType &extractor)
    {
      const std::vector<unsigned int> index_set(
        internal::extract_field_component_indices<dim>(extractor));
      for (unsigned int i = 0; i < index_set.size(); ++i)
        {
          Assert(this->registered_marked_dependent_variables[index_set[i]] ==
                   false,
                 ExcMessage("Overlapping indices for dependent variables."));
          HelperBase<ADNumberTypeCode, ScalarType>::register_dependent_variable(
            index_set[i], internal::get_tensor_entry(funcs, i));
        }
    }



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType>
    template <typename ExtractorType_Row>
    typename internal::VectorFieldValue<
      dim,
      typename VectorFunction<dim, ADNumberTypeCode, ScalarType>::scalar_type,
      ExtractorType_Row>::type
    VectorFunction<dim, ADNumberTypeCode, ScalarType>::extract_value_component(
      const Vector<scalar_type> &values,
      const ExtractorType_Row &  extractor_row)
    {
      // NOTE: The order of components must be consistently defined throughout
      // this class.
      typename internal::VectorFieldValue<dim, scalar_type, ExtractorType_Row>::
        type out;

      // Get indexsets for the subblock from which we wish to extract the
      // gradient values
      const std::vector<unsigned int> row_index_set(
        internal::extract_field_component_indices<dim>(extractor_row));
      Assert(out.n_independent_components == row_index_set.size(),
             ExcMessage("Not all tensor components have been extracted!"));
      for (unsigned int r = 0; r < row_index_set.size(); ++r)
        internal::set_tensor_entry(out, r, values[row_index_set[r]]);

      return out;
    }



    template <int                  dim,
              enum AD::NumberTypes ADNumberTypeCode,
              typename ScalarType>
    template <typename ExtractorType_Row, typename ExtractorType_Col>
    typename internal::VectorFieldJacobian<
      dim,
      typename VectorFunction<dim, ADNumberTypeCode, ScalarType>::scalar_type,
      ExtractorType_Row,
      ExtractorType_Col>::type
    VectorFunction<dim, ADNumberTypeCode, ScalarType>::
      extract_jacobian_component(const FullMatrix<scalar_type> &jacobian,
                                 const ExtractorType_Row &      extractor_row,
                                 const ExtractorType_Col &      extractor_col)
    {
      using InternalJacobian     = internal::VectorFieldJacobian<dim,
                                                             scalar_type,
                                                             ExtractorType_Row,
                                                             ExtractorType_Col>;
      using InternalExtractorRow = internal::Extractor<dim, ExtractorType_Row>;
      using InternalExtractorCol = internal::Extractor<dim, ExtractorType_Col>;
      using JacobianType         = typename InternalJacobian::type;

      // NOTE: The order of components must be consistently defined throughout
      // this class.
      JacobianType out;

      // Get indexsets for the subblocks from which we wish to extract the
      // Hessian values.
      // NOTE: Here we have to do some clever accounting when the
      // one extractor is a symmetric Tensor and the other is not, e.g.
      // <SymmTensor,Vector>. In this scenario the return type is a
      // non-symmetric Tensor<3,dim> but we have to fetch information from a
      // SymmTensor row/column that has too few entries to fill the output
      // tensor. So we must duplicate the relevant entries in the row/column
      // indexset to fetch off-diagonal components that are Otherwise
      // non-existent in a SymmTensor.
      const std::vector<unsigned int> row_index_set(
        internal::extract_field_component_indices<dim>(
          extractor_row, false /*ignore_symmetries*/));
      const std::vector<unsigned int> col_index_set(
        internal::extract_field_component_indices<dim>(
          extractor_col, false /*ignore_symmetries*/));

      for (unsigned int index = 0;
           index < JacobianType::n_independent_components;
           ++index)
        {
          const TableIndices<JacobianType::rank> ti_out =
            JacobianType::unrolled_to_component_indices(index);
          const unsigned int r =
            InternalExtractorRow::local_component(ti_out, 0);
          const unsigned int c =
            InternalExtractorCol::local_component(ti_out,
                                                  InternalExtractorRow::rank);

          internal::set_tensor_entry(
            out, index, jacobian[row_index_set[r]][col_index_set[c]]);
        }

      return out;
    }


  } // namespace AD
} // namespace Differentiation


#  endif // DOXYGEN


DEAL_II_NAMESPACE_CLOSE

#endif // defined(DEAL_II_WITH_ADOLC) || defined(DEAL_II_TRILINOS_WITH_SACADO)

#endif // dealii_differentiation_ad_ad_helpers_h