scratch_data.h

// ---------------------------------------------------------------------
//
// Copyright (C) 2019 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.md at
// the top level directory of deal.II.
//
// ---------------------------------------------------------------------

#ifndef dealii_meshworker_scratch_data_h
#define dealii_meshworker_scratch_data_h

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

#include <deal.II/algorithms/general_data_storage.h>

#include <deal.II/differentiation/ad.h>

#include <deal.II/fe/fe_values.h>
#include <deal.II/fe/mapping_q1.h>

#include <boost/any.hpp>

#include <algorithm>
#include <map>
#include <vector>

DEAL_II_NAMESPACE_OPEN

namespace MeshWorker
{
  namespace internal
  {
    template <int dim, int spacedim, typename NumberType, typename Extractor>
    using OutputType = typename FEValuesViews::View<dim, spacedim, Extractor>::
      template OutputType<NumberType>;
  } // namespace internal

  template <int dim, int spacedim = dim>
  class ScratchData
  {
  public:
    ScratchData(
      const Mapping<dim, spacedim> &      mapping,
      const FiniteElement<dim, spacedim> &fe,
      const Quadrature<dim> &             quadrature,
      const UpdateFlags &                 update_flags,
      const Quadrature<dim - 1> &face_quadrature   = Quadrature<dim - 1>(),
      const UpdateFlags &        face_update_flags = update_default);

    ScratchData(
      const Mapping<dim, spacedim> &      mapping,
      const FiniteElement<dim, spacedim> &fe,
      const Quadrature<dim> &             quadrature,
      const UpdateFlags &                 update_flags,
      const UpdateFlags &                 neighbor_update_flags,
      const Quadrature<dim - 1> &face_quadrature   = Quadrature<dim - 1>(),
      const UpdateFlags &        face_update_flags = update_default,
      const UpdateFlags &        neighbor_face_update_flags = update_default);

    ScratchData(
      const FiniteElement<dim, spacedim> &fe,
      const Quadrature<dim> &             quadrature,
      const UpdateFlags &                 update_flags,
      const Quadrature<dim - 1> &face_quadrature   = Quadrature<dim - 1>(),
      const UpdateFlags &        face_update_flags = update_default);

    ScratchData(
      const FiniteElement<dim, spacedim> &fe,
      const Quadrature<dim> &             quadrature,
      const UpdateFlags &                 update_flags,
      const UpdateFlags &                 neighbor_update_flags,
      const Quadrature<dim - 1> &face_quadrature   = Quadrature<dim - 1>(),
      const UpdateFlags &        face_update_flags = update_default,
      const UpdateFlags &        neighbor_face_update_flags = update_default);

    ScratchData(const ScratchData<dim, spacedim> &scratch);
 // CurrentCellMethods

    const FEValues<dim, spacedim> &
    reinit(
      const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell);

    const FEFaceValues<dim, spacedim> &
    reinit(const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
           const unsigned int face_no);

    const FEFaceValuesBase<dim, spacedim> &
    reinit(const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
           const unsigned int face_no,
           const unsigned int subface_no);

    const FEValuesBase<dim, spacedim> &
    get_current_fe_values() const;

    const std::vector<Point<spacedim>> &
    get_quadrature_points() const;

    const std::vector<double> &
    get_JxW_values() const;

    const std::vector<Tensor<1, spacedim>> &
    get_normal_vectors() const;

    const std::vector<types::global_dof_index> &
    get_local_dof_indices() const;
 // CurrentCellMethods
 // NeighborCellMethods

    const FEValues<dim, spacedim> &
    reinit_neighbor(
      const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell);

    const FEFaceValues<dim, spacedim> &
    reinit_neighbor(
      const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
      const unsigned int                                              face_no);

    const FEFaceValuesBase<dim, spacedim> &
    reinit_neighbor(
      const typename DoFHandler<dim, spacedim>::active_cell_iterator &cell,
      const unsigned int                                              face_no,
      const unsigned int subface_no);

    const FEValuesBase<dim, spacedim> &
    get_current_neighbor_fe_values() const;

    const std::vector<double> &
    get_neighbor_JxW_values() const;

    const std::vector<Tensor<1, spacedim>> &
    get_neighbor_normal_vectors();

    const std::vector<types::global_dof_index> &
    get_neighbor_dof_indices() const;
 // NeighborCellMethods

    GeneralDataStorage &
    get_general_data_storage();
 // CurrentCellEvaluation

    template <typename VectorType, typename Number = double>
    void
    extract_local_dof_values(const std::string &global_vector_name,
                             const VectorType & input_vector,
                             const Number       dummy = Number(0));

    template <typename Number = double>
    const std::vector<Number> &
    get_local_dof_values(const std::string &global_vector_name,
                         Number             dummy = Number(0)) const;

    template <typename Extractor, typename Number = double>
    const std::vector<
      typename internal::OutputType<dim, spacedim, Number, Extractor>::
        value_type> &
    get_values(const std::string &global_vector_name,
               const Extractor &  variable,
               const Number       dummy = Number(0));

    template <typename Extractor, typename Number = double>
    const std::vector<
      typename internal::OutputType<dim, spacedim, Number, Extractor>::
        gradient_type> &
    get_gradients(const std::string &global_vector_name,
                  const Extractor &  variable,
                  const Number       dummy = Number(0));

    template <typename Extractor, typename Number = double>
    const std::vector<
      typename internal::OutputType<dim, spacedim, Number, Extractor>::
        symmetric_gradient_type> &
    get_symmetric_gradients(const std::string &global_vector_name,
                            const Extractor &  variable,
                            const Number       dummy = Number(0));

    template <typename Extractor, typename Number = double>
    const std::vector<
      typename internal::OutputType<dim, spacedim, Number, Extractor>::
        divergence_type> &
    get_divergences(const std::string &global_vector_name,
                    const Extractor &  variable,
                    const Number       dummy = Number(0));

    template <typename Extractor, typename Number = double>
    const std::vector<typename internal::
                        OutputType<dim, spacedim, Number, Extractor>::curl_type>
      &
      get_curls(const std::string &global_vector_name,
                const Extractor &  variable,
                const Number       dummy = Number(0));

    template <typename Extractor, typename Number = double>
    const std::vector<
      typename internal::OutputType<dim, spacedim, Number, Extractor>::
        hessian_type> &
    get_hessians(const std::string &global_vector_name,
                 const Extractor &  variable,
                 const Number       dummy = Number(0));

    template <typename Extractor, typename Number = double>
    const std::vector<
      typename internal::OutputType<dim, spacedim, Number, Extractor>::
        third_derivative_type> &
    get_third_derivatives(const std::string &global_vector_name,
                          const Extractor &  variable,
                          const Number       dummy = Number(0));
 // CurrentCellEvaluation

    const Mapping<dim, spacedim> &
    get_mapping() const;

  private:
    template <typename Extractor, typename Number = double>
    std::string
    get_unique_name(const std::string &global_vector_name,
                    const Extractor &  variable,
                    const std::string &object_type,
                    const unsigned int size,
                    const Number &     exemplar_number) const;

    template <typename Number = double>
    std::string
    get_unique_dofs_name(const std::string &global_vector_name,
                         const unsigned int size,
                         const Number &     exemplar_number) const;

    SmartPointer<const Mapping<dim, spacedim>> mapping;

    SmartPointer<const FiniteElement<dim, spacedim>> fe;

    Quadrature<dim> cell_quadrature;

    Quadrature<dim - 1> face_quadrature;

    UpdateFlags cell_update_flags;

    UpdateFlags neighbor_cell_update_flags;

    UpdateFlags face_update_flags;

    UpdateFlags neighbor_face_update_flags;

    std::unique_ptr<FEValues<dim, spacedim>> fe_values;

    std::unique_ptr<FEFaceValues<dim, spacedim>> fe_face_values;

    std::unique_ptr<FESubfaceValues<dim, spacedim>> fe_subface_values;

    std::unique_ptr<FEValues<dim, spacedim>> neighbor_fe_values;

    std::unique_ptr<FEFaceValues<dim, spacedim>> neighbor_fe_face_values;

    std::unique_ptr<FESubfaceValues<dim, spacedim>> neighbor_fe_subface_values;

    std::vector<types::global_dof_index> local_dof_indices;

    std::vector<types::global_dof_index> neighbor_dof_indices;

    GeneralDataStorage user_data_storage;

    GeneralDataStorage internal_data_storage;

    SmartPointer<FEValuesBase<dim, spacedim>> current_fe_values;

    SmartPointer<FEValuesBase<dim, spacedim>> current_neighbor_fe_values;
  };

#ifndef DOXYGEN
  template <int dim, int spacedim>
  template <typename Extractor, typename Number>
  std::string
  ScratchData<dim, spacedim>::get_unique_name(
    const std::string &global_vector_name,
    const Extractor &  variable,
    const std::string &object_type,
    const unsigned int size,
    const Number &     exemplar_number) const
  {
    return global_vector_name + "_" + variable.get_name() + "_" + object_type +
           "_" + Utilities::int_to_string(size) + "_" +
           Utilities::type_to_string(exemplar_number);
  }



  template <int dim, int spacedim>
  template <typename Number>
  std::string
  ScratchData<dim, spacedim>::get_unique_dofs_name(
    const std::string &global_vector_name,
    const unsigned int size,
    const Number &     exemplar_number) const
  {
    return global_vector_name + "_independent_local_dofs_" +
           Utilities::int_to_string(size) + "_" +
           Utilities::type_to_string(exemplar_number);
  }



  template <int dim, int spacedim>
  template <typename VectorType, typename Number>
  void
  ScratchData<dim, spacedim>::extract_local_dof_values(
    const std::string &global_vector_name,
    const VectorType & input_vector,
    const Number       dummy)
  {
    const unsigned int n_dofs = get_current_fe_values().get_fe().dofs_per_cell;

    const std::string name =
      get_unique_dofs_name(global_vector_name, n_dofs, dummy);

    auto &independent_local_dofs =
      internal_data_storage
        .template get_or_add_object_with_name<std::vector<Number>>(name,
                                                                   n_dofs);

    AssertDimension(independent_local_dofs.size(), n_dofs);

    if (Differentiation::AD::is_tapeless_ad_number<Number>::value == true)
      for (unsigned int i = 0; i < n_dofs; ++i)
        Differentiation::AD::internal::Marking<Number>::independent_variable(
          input_vector(local_dof_indices[i]),
          i,
          n_dofs,
          independent_local_dofs[i]);
    else
      for (unsigned int i = 0; i < n_dofs; ++i)
        independent_local_dofs[i] = input_vector(local_dof_indices[i]);
  }



  template <int dim, int spacedim>
  template <typename Number>
  const std::vector<Number> &
  ScratchData<dim, spacedim>::get_local_dof_values(
    const std::string &global_vector_name,
    Number             dummy) const
  {
    const unsigned int n_dofs = get_current_fe_values().get_fe().dofs_per_cell;

    const std::string dofs_name =
      get_unique_dofs_name(global_vector_name, n_dofs, dummy);

    Assert(
      internal_data_storage.stores_object_with_name(dofs_name),
      ExcMessage(
        "You did not call yet extract_local_dof_values with the right types!"));

    return internal_data_storage
      .template get_object_with_name<std::vector<Number>>(dofs_name);
  }



  template <int dim, int spacedim>
  template <typename Extractor, typename Number>
  const std::vector<
    typename internal::OutputType<dim, spacedim, Number, Extractor>::value_type>
    &
    ScratchData<dim, spacedim>::get_values(
      const std::string &global_vector_name,
      const Extractor &  variable,
      const Number       dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_local_dof_values(global_vector_name, dummy);

    const FEValuesBase<dim, spacedim> &fev = get_current_fe_values();

    const unsigned int n_q_points = fev.n_quadrature_points;

    const std::string name = get_unique_name(
      global_vector_name, variable, "_values_q", n_q_points, dummy);

    // Now build the return type
    using RetType =
      std::vector<typename internal::
                    OutputType<dim, spacedim, Number, Extractor>::value_type>;

    RetType &ret =
      internal_data_storage.template get_or_add_object_with_name<RetType>(
        name, n_q_points);

    AssertDimension(ret.size(), n_q_points);

    fev[variable].get_function_values_from_local_dof_values(
      independent_local_dofs, ret);
    return ret;
  }



  template <int dim, int spacedim>
  template <typename Extractor, typename Number>
  const std::vector<
    typename internal::OutputType<dim, spacedim, Number, Extractor>::
      gradient_type> &
  ScratchData<dim, spacedim>::get_gradients(
    const std::string &global_vector_name,
    const Extractor &  variable,
    const Number       dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_local_dof_values(global_vector_name, dummy);

    const FEValuesBase<dim, spacedim> &fev = get_current_fe_values();

    const unsigned int n_q_points = fev.n_quadrature_points;

    const std::string name = get_unique_name(
      global_vector_name, variable, "_gradients_q", n_q_points, dummy);

    // Now build the return type
    using RetType = std::vector<
      typename internal::OutputType<dim, spacedim, Number, Extractor>::
        gradient_type>;

    RetType &ret =
      internal_data_storage.template get_or_add_object_with_name<RetType>(
        name, n_q_points);

    AssertDimension(ret.size(), n_q_points);

    fev[variable].get_function_gradients_from_local_dof_values(
      independent_local_dofs, ret);
    return ret;
  }



  template <int dim, int spacedim>
  template <typename Extractor, typename Number>
  const std::vector<
    typename internal::OutputType<dim, spacedim, Number, Extractor>::
      hessian_type> &
  ScratchData<dim, spacedim>::get_hessians(
    const std::string &global_vector_name,
    const Extractor &  variable,
    const Number       dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_local_dof_values(global_vector_name, dummy);

    const FEValuesBase<dim, spacedim> &fev = get_current_fe_values();

    const unsigned int n_q_points = fev.n_quadrature_points;

    const std::string name = get_unique_name(
      global_vector_name, variable, "_hessians_q", n_q_points, dummy);

    // Now build the return type
    using RetType =
      std::vector<typename internal::
                    OutputType<dim, spacedim, Number, Extractor>::hessian_type>;

    RetType &ret =
      internal_data_storage.template get_or_add_object_with_name<RetType>(
        name, n_q_points);


    AssertDimension(ret.size(), n_q_points);

    fev[variable].get_function_hessians_from_local_dof_values(
      independent_local_dofs, ret);
    return ret;
  }



  template <int dim, int spacedim>
  template <typename Extractor, typename Number>
  const std::vector<
    typename internal::OutputType<dim, spacedim, Number, Extractor>::
      third_derivative_type> &
  ScratchData<dim, spacedim>::get_third_derivatives(
    const std::string &global_vector_name,
    const Extractor &  variable,
    const Number       dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_local_dof_values(global_vector_name, dummy);

    const FEValuesBase<dim, spacedim> &fev = get_current_fe_values();

    const unsigned int n_q_points = fev.n_quadrature_points;

    const std::string name = get_unique_name(
      global_vector_name, variable, "_third_derivatives_q", n_q_points, dummy);

    // Now build the return type
    using RetType = std::vector<
      typename internal::OutputType<dim, spacedim, Number, Extractor>::
        third_derivative_type>;

    RetType &ret =
      internal_data_storage.template get_or_add_object_with_name<RetType>(
        name, n_q_points);


    AssertDimension(ret.size(), n_q_points);

    fev[variable].get_function_third_derivatives_from_local_dof_values(
      independent_local_dofs, ret);
    return ret;
  }



  template <int dim, int spacedim>
  template <typename Extractor, typename Number>
  const std::vector<
    typename internal::OutputType<dim, spacedim, Number, Extractor>::
      symmetric_gradient_type> &
  ScratchData<dim, spacedim>::get_symmetric_gradients(
    const std::string &global_vector_name,
    const Extractor &  variable,
    const Number       dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_local_dof_values(global_vector_name, dummy);

    const FEValuesBase<dim, spacedim> &fev = get_current_fe_values();

    const unsigned int n_q_points = fev.n_quadrature_points;

    const std::string name = get_unique_name(
      global_vector_name, variable, "_symmetric_gradient_q", n_q_points, dummy);


    // Now build the return type
    using RetType = std::vector<
      typename internal::OutputType<dim, spacedim, Number, Extractor>::
        symmetric_gradient_type>;

    RetType &ret =
      internal_data_storage.template get_or_add_object_with_name<RetType>(
        name, n_q_points);


    AssertDimension(ret.size(), n_q_points);

    fev[variable].get_function_symmetric_gradients_from_local_dof_values(
      independent_local_dofs, ret);
    return ret;
  }


  template <int dim, int spacedim>
  template <typename Extractor, typename Number>
  const std::vector<
    typename internal::OutputType<dim, spacedim, Number, Extractor>::
      divergence_type> &
  ScratchData<dim, spacedim>::get_divergences(
    const std::string &global_vector_name,
    const Extractor &  variable,
    const Number       dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_local_dof_values(global_vector_name, dummy);

    const FEValuesBase<dim, spacedim> &fev = get_current_fe_values();

    const unsigned int n_q_points = fev.n_quadrature_points;

    const std::string name = get_unique_name(
      global_vector_name, variable, "_divergence_q", n_q_points, dummy);

    // Now build the return type
    using RetType = std::vector<
      typename internal::OutputType<dim, spacedim, Number, Extractor>::
        divergence_type>;

    RetType &ret =
      internal_data_storage.template get_or_add_object_with_name<RetType>(
        name, n_q_points);


    AssertDimension(ret.size(), n_q_points);

    fev[variable].get_function_divergences_from_local_dof_values(
      independent_local_dofs, ret);
    return ret;
  }



  template <int dim, int spacedim>
  template <typename Extractor, typename Number>
  const std::vector<
    typename internal::OutputType<dim, spacedim, Number, Extractor>::curl_type>
    &
    ScratchData<dim, spacedim>::get_curls(const std::string &global_vector_name,
                                          const Extractor &  variable,
                                          const Number       dummy)
  {
    const std::vector<Number> &independent_local_dofs =
      get_local_dof_values(global_vector_name, dummy);

    const FEValuesBase<dim, spacedim> &fev = get_current_fe_values();

    const unsigned int n_q_points = fev.n_quadrature_points;

    const std::string name = get_unique_name(
      global_vector_name, variable, "_curl_q", n_q_points, dummy);

    // Now build the return type
    using RetType =
      std::vector<typename internal::
                    OutputType<dim, spacedim, Number, Extractor>::curl_type>;

    RetType &ret =
      internal_data_storage.template get_or_add_object_with_name<RetType>(
        name, n_q_points);

    AssertDimension(ret.size(), n_q_points);

    fev[variable].get_function_curls_from_local_dof_values(
      independent_local_dofs, ret);
    return ret;
  }

#endif

} // namespace MeshWorker

DEAL_II_NAMESPACE_CLOSE

#endif