fe_coupling_values.h

Go to the documentation of this file.

// ---------------------------------------------------------------------
//
// Copyright (C) 2023 - 2024 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_fe_coupling_values_h
#define dealii_fe_coupling_values_h
 
#include <deal.II/base/config.h>
 
#include <deal.II/algorithms/general_data_storage.h>
 
#include <deal.II/base/enable_observer_pointer.h>
#include <deal.II/base/std_cxx20/iota_view.h>
#include <deal.II/base/thread_local_storage.h>
#include <deal.II/base/utilities.h>
 
#include <deal.II/fe/fe_values_views.h>
 
DEAL_II_NAMESPACE_OPEN
 
// Forward declaration
#ifndef DOXYGEN
template <int dim, int spacedim>
class FEValuesBase;
#endif
 
namespace FEValuesViews
{
  struct RenumberingData
  {
    RenumberingData(
      const unsigned int n_inner_dofs = numbers::invalid_unsigned_int,
      const unsigned int n_inner_quadrature_points =
        numbers::invalid_unsigned_int,
      const std::vector<unsigned int> &dof_renumbering        = {},
      const std::vector<unsigned int> &quadrature_renumbering = {});
 
    RenumberingData(const RenumberingData &other) = delete;
 
    const unsigned int n_inner_dofs;
 
    const unsigned int n_dofs;
 
    const unsigned int n_inner_quadrature_points;
 
    const unsigned int n_quadrature_points;
 
    const std::vector<unsigned int> dof_renumbering;
 
    const std::vector<unsigned int> quadrature_renumbering;
 
    mutable Threads::ThreadLocalStorage<GeneralDataStorage> data_storage;
  };
 
  template <typename ViewType>
  class RenumberedView
  {
  public:
    using value_type = typename ViewType::value_type;
 
    using gradient_type = typename ViewType::gradient_type;
 
    template <typename Number>
    using solution_value_type =
      typename ViewType::template solution_value_type<Number>;
 
    template <typename Number>
    using solution_gradient_type =
      typename ViewType::template solution_gradient_type<Number>;
 
    RenumberedView(const ViewType &view, const RenumberingData &data);
 
    value_type
    value(const unsigned int shape_function, const unsigned int q_point) const;
 
    gradient_type
    gradient(const unsigned int shape_function,
             const unsigned int q_point) const;
 
    template <typename Number>
    void
    get_function_values(const ReadVector<Number>                 &fe_function,
                        std::vector<solution_value_type<Number>> &values) const;
 
    template <class InputVector>
    void
    get_function_values_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_value_type<typename InputVector::value_type>>
        &values) const;
 
    template <typename Number>
    void
    get_function_gradients(
      const ReadVector<Number>                    &fe_function,
      std::vector<solution_gradient_type<Number>> &gradients) const;
 
    template <class InputVector>
    void
    get_function_gradients_from_local_dof_values(
      const InputVector &dof_values,
      std::vector<solution_gradient_type<typename InputVector::value_type>>
        &gradients) const;
 
  private:
    const RenumberingData &data;
 
 
    template <typename Number>
    std::string
    get_unique_container_name(const std::string &prefix,
                              const unsigned int size,
                              const Number      &exemplar_number) const;
 
    template <typename InputVector>
    const InputVector &
    outer_to_inner_dofs(const InputVector &outer_vector) const;
 
    template <typename ValueType>
    std::vector<ValueType> &
    outer_to_inner_values(std::vector<ValueType> &outer_values) const;
 
    template <typename ValueType>
    void
    inner_to_outer_values(const std::vector<ValueType> &inner_values,
                          std::vector<ValueType>       &outer_values) const;
 
    const ViewType &view;
  };
} // namespace FEValuesViews
 
 
enum class QuadratureCouplingType
{
  tensor_product,
 
  unrolled,
 
  matching,
 
  reorder,
 
  overlapping,
};
 
enum class DoFCouplingType
{
  independent,
 
  contiguous,
};
 
template <int dim1, int dim2 = dim1, int spacedim = dim1>
class FECouplingValues
{
public:
  FECouplingValues();
 
  FECouplingValues(
    const FEValuesBase<dim1, spacedim> &fe_values_1,
    const FEValuesBase<dim2, spacedim> &fe_values_2,
    const DoFCouplingType &dof_coupling_type = DoFCouplingType::independent,
    const QuadratureCouplingType &quadrature_coupling_type =
      QuadratureCouplingType::tensor_product);
 
  void
  reinit(
    const FEValuesBase<dim1, spacedim> &fe_values_1,
    const FEValuesBase<dim2, spacedim> &fe_values_2,
    const DoFCouplingType &dof_coupling_type = DoFCouplingType::independent,
    const QuadratureCouplingType &quadrature_coupling_type =
      QuadratureCouplingType::tensor_product);
 
 
  template <typename Extractor>
  FEValuesExtractors::FirstCoupling<Extractor>
  get_first_extractor(const Extractor &extractor) const;
 
  template <typename Extractor>
  FEValuesExtractors::SecondCoupling<Extractor>
  get_second_extractor(const Extractor &extractor) const;
 
  double
  JxW(const unsigned int quadrature_point) const;
 
  std::pair<Point<spacedim>, Point<spacedim>>
  quadrature_point(const unsigned int quadrature_point) const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  quadrature_point_indices() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  first_dof_indices() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  second_dof_indices() const;
 
  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
  coupling_dof_indices() const;
 
  std::vector<types::global_dof_index>
  get_coupling_dof_indices(
    const std::vector<types::global_dof_index> &dof_indices_1,
    const std::vector<types::global_dof_index> &dof_indices_2,
    const types::global_dof_index               dofs_offset_1 = 0,
    const types::global_dof_index               dofs_offset_2 = 0) const;
 
  std::pair<unsigned int, unsigned int>
  coupling_dof_to_dof_indices(const unsigned int coupling_dof_index) const;
 
  std::pair<unsigned int, unsigned int>
  coupling_quadrature_to_quadrature_indices(
    const unsigned int quadrature_point) const;
 
  template <typename Extractor>
  const FEValuesViews::RenumberedView<
    FEValuesViews::View<dim1, spacedim, Extractor>>
  operator[](
    const FEValuesExtractors::FirstCoupling<Extractor> &extractor) const;
 
  template <typename Extractor>
  const FEValuesViews::RenumberedView<
    FEValuesViews::View<dim2, spacedim, Extractor>>
  operator[](
    const FEValuesExtractors::SecondCoupling<Extractor> &extractor) const;
 
  unsigned int
  n_coupling_dofs() const;
 
  unsigned int
  n_first_dofs() const;
 
  unsigned int
  n_second_dofs() const;
 
  unsigned int
  n_quadrature_points() const;
 
private:
  DoFCouplingType dof_coupling_type;
 
  QuadratureCouplingType quadrature_coupling_type;
 
  ObserverPointer<const FEValuesBase<dim1, spacedim>> first_fe_values;
 
  ObserverPointer<const FEValuesBase<dim2, spacedim>> second_fe_values;
 
  std::unique_ptr<const FEValuesViews::RenumberingData> first_renumbering_data;
 
  std::unique_ptr<const FEValuesViews::RenumberingData> second_renumbering_data;
 
  unsigned int n_quadrature_points_;
 
  unsigned int n_coupling_dofs_;
};
 
 
#ifndef DOXYGEN
 
 
/*------------------------ Inline functions: namespace FEValuesViews --------*/
 
namespace FEValuesViews
{
  RenumberingData::RenumberingData(
    const unsigned int               n_inner_dofs,
    const unsigned int               n_inner_quadrature_points,
    const std::vector<unsigned int> &dof_renumbering,
    const std::vector<unsigned int> &quadrature_renumbering)
    : n_inner_dofs(n_inner_dofs)
    , n_dofs(dof_renumbering.empty() ? n_inner_dofs : dof_renumbering.size())
    , n_inner_quadrature_points(n_inner_quadrature_points)
    , n_quadrature_points(quadrature_renumbering.empty() ?
                            n_inner_quadrature_points :
                            quadrature_renumbering.size())
    , dof_renumbering(dof_renumbering)
    , quadrature_renumbering(quadrature_renumbering)
  {
// Check that the renumbering vectors are valid.
#  ifdef DEBUG
    // While for dofs we admit invalid values, this is not the case for
    // quadrature points.
    for (const auto i : dof_renumbering)
      Assert(i < n_inner_dofs || i == numbers::invalid_unsigned_int,
             ExcIndexRange(i, 0, n_inner_dofs));
 
    for (const auto q : quadrature_renumbering)
      AssertIndexRange(q, n_inner_quadrature_points);
#  endif
  }
 
 
 
  template <typename ViewType>
  RenumberedView<ViewType>::RenumberedView(const ViewType        &view,
                                           const RenumberingData &data)
    : view(view)
    , data(data)
  {}
 
 
 
  template <typename ViewType>
  template <typename Number>
  inline std::string
  RenumberedView<ViewType>::get_unique_container_name(
    const std::string &prefix,
    const unsigned int size,
    const Number      &exemplar_number) const
  {
    return prefix + "_" + Utilities::int_to_string(size) + "_" +
           Utilities::type_to_string(exemplar_number);
  }
 
 
 
  template <typename ViewType>
  typename RenumberedView<ViewType>::value_type
  RenumberedView<ViewType>::value(const unsigned int shape_function,
                                  const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, data.n_dofs);
    AssertIndexRange(q_point, data.n_quadrature_points);
 
    const auto inner_shape_function = data.dof_renumbering.empty() ?
                                        shape_function :
                                        data.dof_renumbering[shape_function];
    const auto inner_q_point        = data.quadrature_renumbering.empty() ?
                                        q_point :
                                        data.quadrature_renumbering[q_point];
    if (inner_shape_function == numbers::invalid_unsigned_int)
      return value_type(0);
    else
      {
        AssertIndexRange(inner_shape_function, data.n_inner_dofs);
        AssertIndexRange(inner_q_point, data.n_inner_quadrature_points);
        return view.value(inner_shape_function, inner_q_point);
      }
  }
 
 
 
  template <typename ViewType>
  typename RenumberedView<ViewType>::gradient_type
  RenumberedView<ViewType>::gradient(const unsigned int shape_function,
                                     const unsigned int q_point) const
  {
    AssertIndexRange(shape_function, data.n_dofs);
    AssertIndexRange(q_point, data.n_quadrature_points);
 
    const auto inner_shape_function = data.dof_renumbering.empty() ?
                                        shape_function :
                                        data.dof_renumbering[shape_function];
    const auto inner_q_point        = data.quadrature_renumbering.empty() ?
                                        q_point :
                                        data.quadrature_renumbering[q_point];
    if (inner_shape_function == numbers::invalid_unsigned_int)
      return gradient_type();
    else
      return view.gradient(inner_shape_function, inner_q_point);
  }
 
 
 
  template <typename ViewType>
  template <typename ValueType>
  std::vector<ValueType> &
  RenumberedView<ViewType>::outer_to_inner_values(
    std::vector<ValueType> &outer_values) const
  {
    AssertDimension(outer_values.size(), data.n_quadrature_points);
    if (data.quadrature_renumbering.empty())
      {
        return outer_values;
      }
    else
      {
        const auto name =
          get_unique_container_name("RenumberedView::outer_to_inner_values",
                                    data.n_inner_quadrature_points,
                                    outer_values[0]);
        auto &inner_values =
          data.data_storage.get()
            .template get_or_add_object_with_name<std::vector<ValueType>>(
              name, data.n_inner_quadrature_points);
        return inner_values;
      }
  }
 
 
 
  template <typename ViewType>
  template <typename VectorType>
  const VectorType &
  RenumberedView<ViewType>::outer_to_inner_dofs(
    const VectorType &outer_dofs) const
  {
    AssertDimension(outer_dofs.size(), data.n_dofs);
    if (data.dof_renumbering.empty())
      {
        return outer_dofs;
      }
    else
      {
        const auto name =
          get_unique_container_name("RenumberedView::outer_to_inner_dofs",
                                    data.n_inner_dofs,
                                    outer_dofs[0]);
 
        auto &inner_dofs = data.data_storage.get()
                             .template get_or_add_object_with_name<VectorType>(
                               name, data.n_inner_dofs);
        for (unsigned int i = 0; i < data.n_dofs; ++i)
          {
            const auto inner_i = data.dof_renumbering[i];
            if (inner_i != numbers::invalid_unsigned_int)
              {
                AssertIndexRange(inner_i, data.n_inner_dofs);
                inner_dofs[inner_i] = outer_dofs[i];
              }
          }
        return inner_dofs;
      }
  }
 
 
 
  template <typename ViewType>
  template <typename ValueType>
  void
  RenumberedView<ViewType>::inner_to_outer_values(
    const std::vector<ValueType> &inner_values,
    std::vector<ValueType>       &outer_values) const
  {
    AssertDimension(outer_values.size(), data.n_quadrature_points);
    AssertDimension(inner_values.size(), data.n_inner_quadrature_points);
    if (data.quadrature_renumbering.empty())
      {
        Assert(&inner_values == &outer_values, ExcInternalError());
        return;
      }
    for (unsigned int i = 0; i < data.quadrature_renumbering.size(); ++i)
      {
        outer_values[i] = inner_values[data.quadrature_renumbering[i]];
      }
  }
 
 
 
  template <typename ViewType>
  template <typename Number>
  void
  RenumberedView<ViewType>::get_function_values(
    const ReadVector<Number>                 &fe_function,
    std::vector<solution_value_type<Number>> &values) const
  {
    auto &inner_values = outer_to_inner_values(values);
    view.get_function_values(fe_function, inner_values);
    inner_to_outer_values(inner_values, values);
  }
 
 
 
  template <typename ViewType>
  template <typename InputVector>
  void
  RenumberedView<ViewType>::get_function_values_from_local_dof_values(
    const InputVector &dof_values,
    std::vector<solution_value_type<typename InputVector::value_type>> &values)
    const
  {
    const auto &inner_dof_values = outer_to_inner_dofs(dof_values);
    auto       &inner_values     = outer_to_inner_values(values);
 
    view.get_function_values_from_local_dof_values(inner_dof_values,
                                                   inner_values);
    inner_to_outer_values(inner_values, values);
  }
 
 
  template <typename ViewType>
  template <typename Number>
  void
  RenumberedView<ViewType>::get_function_gradients(
    const ReadVector<Number>                    &fe_function,
    std::vector<solution_gradient_type<Number>> &gradients) const
  {
    auto &inner_gradients = outer_to_inner_values(gradients);
    view.get_function_gradients(fe_function, inner_gradients);
    inner_to_outer_values(inner_gradients, gradients);
  }
 
 
 
  template <typename ViewType>
  template <typename InputVector>
  void
  RenumberedView<ViewType>::get_function_gradients_from_local_dof_values(
    const InputVector &dof_values,
    std::vector<solution_gradient_type<typename InputVector::value_type>>
      &gradients) const
  {
    const auto &inner_dof_values = outer_to_inner_dofs(dof_values);
    auto       &inner_gradients  = outer_to_inner_values(gradients);
 
    view.get_function_gradients_from_local_dof_values(inner_dof_values,
                                                      inner_gradients);
    inner_to_outer_values(inner_gradients, gradients);
  }
} // namespace FEValuesViews
 
/*-------------- Inline functions FECouplingValues ---------------------*/
 
template <int dim1, int dim2, int spacedim>
std::vector<types::global_dof_index>
FECouplingValues<dim1, dim2, spacedim>::get_coupling_dof_indices(
  const std::vector<types::global_dof_index> &dof_indices_1,
  const std::vector<types::global_dof_index> &dof_indices_2,
  const types::global_dof_index               dofs_offset_1,
  const types::global_dof_index               dofs_offset_2) const
{
  AssertThrow(
    n_coupling_dofs_ != numbers::invalid_unsigned_int,
    ExcMessage(
      "Dofs are independent. You cannot ask for coupling dof indices."));
  AssertDimension(dof_indices_1.size(), first_fe_values->dofs_per_cell);
  AssertDimension(dof_indices_2.size(), second_fe_values->dofs_per_cell);
 
  std::vector<types::global_dof_index> coupling_dof_indices(
    dof_indices_1.size() + dof_indices_2.size());
  unsigned int idx = 0;
  for (const auto &i : dof_indices_1)
    coupling_dof_indices[idx++] = i + dofs_offset_1;
  for (const auto &i : dof_indices_2)
    coupling_dof_indices[idx++] = i + dofs_offset_2;
  return coupling_dof_indices;
}
 
 
 
template <int dim1, int dim2, int spacedim>
FECouplingValues<dim1, dim2, spacedim>::FECouplingValues()
  : first_fe_values(nullptr)
  , second_fe_values(nullptr)
  , quadrature_coupling_type(QuadratureCouplingType::unrolled)
  , n_quadrature_points_(0)
  , n_coupling_dofs_(numbers::invalid_unsigned_int)
{}
 
 
 
template <int dim1, int dim2, int spacedim>
FECouplingValues<dim1, dim2, spacedim>::FECouplingValues(
  const FEValuesBase<dim1, spacedim> &fe_values_1,
  const FEValuesBase<dim2, spacedim> &fe_values_2,
  const DoFCouplingType              &dof_coupling_type,
  const QuadratureCouplingType       &quadrature_coupling_type)
  : FECouplingValues() // delegate to other constructor
{
  reinit(fe_values_1, fe_values_2, dof_coupling_type, quadrature_coupling_type);
}
 
 
 
template <int dim1, int dim2, int spacedim>
void
FECouplingValues<dim1, dim2, spacedim>::reinit(
  const FEValuesBase<dim1, spacedim> &fe_values_1,
  const FEValuesBase<dim2, spacedim> &fe_values_2,
  const DoFCouplingType              &dof_coupling_type,
  const QuadratureCouplingType       &quadrature_coupling_type)
{
  first_fe_values                = &fe_values_1;
  second_fe_values               = &fe_values_2;
  this->dof_coupling_type        = dof_coupling_type;
  this->quadrature_coupling_type = quadrature_coupling_type;
 
  // Gather information about the inner objects
  unsigned int first_n_inner_dofs  = fe_values_1.dofs_per_cell;
  unsigned int second_n_inner_dofs = fe_values_2.dofs_per_cell;
 
  unsigned int first_n_inner_quadrature_points =
    fe_values_1.n_quadrature_points;
  unsigned int second_n_inner_quadrature_points =
    fe_values_2.n_quadrature_points;
 
  // Initialize counters and renumbering vectors to zero
  unsigned int first_n_dofs  = 0;
  unsigned int second_n_dofs = 0;
 
  unsigned int first_n_quadrature_points  = 0;
  unsigned int second_n_quadrature_points = 0;
 
  std::vector<unsigned int> first_dofs_map;
  std::vector<unsigned int> second_dofs_map;
 
  std::vector<unsigned int> first_quad_map;
  std::vector<unsigned int> second_quad_map;
 
  // Local relative tolerance for comparing quadrature points
  const double tol = std::min(fe_values_1.get_cell()->diameter(),
                              fe_values_2.get_cell()->diameter()) *
                     1e-6;
 
  // Now fill the renumbering vectors
  switch (dof_coupling_type)
    {
      case DoFCouplingType::independent:
        {
          n_coupling_dofs_ = numbers::invalid_unsigned_int;
          first_dofs_map.clear();
          second_dofs_map.clear();
          first_n_dofs  = first_fe_values->dofs_per_cell;
          second_n_dofs = second_fe_values->dofs_per_cell;
          break;
        }
      case DoFCouplingType::contiguous:
        {
          n_coupling_dofs_ =
            fe_values_1.dofs_per_cell + fe_values_2.dofs_per_cell;
 
          first_n_dofs  = n_coupling_dofs_;
          second_n_dofs = n_coupling_dofs_;
 
          first_dofs_map.resize(n_coupling_dofs_);
          second_dofs_map.resize(n_coupling_dofs_);
 
          unsigned int idx = 0;
          for (const unsigned int &i : fe_values_1.dof_indices())
            {
              first_dofs_map[idx]    = i;
              second_dofs_map[idx++] = numbers::invalid_unsigned_int;
            }
          for (const unsigned int &i : fe_values_2.dof_indices())
            {
              first_dofs_map[idx]    = numbers::invalid_unsigned_int;
              second_dofs_map[idx++] = i;
            }
          // Make sure we have the right number of dofs
          AssertDimension(idx, n_coupling_dofs_);
          break;
        }
      default:
        AssertThrow(false, ExcNotImplemented());
    }
  // Compute the quadrature maps
  switch (quadrature_coupling_type)
    {
      case QuadratureCouplingType::tensor_product:
        {
          const auto &quadrature_points_1 = fe_values_1.get_quadrature_points();
          const auto &quadrature_points_2 = fe_values_2.get_quadrature_points();
 
          n_quadrature_points_ =
            quadrature_points_1.size() * quadrature_points_2.size();
 
          first_n_quadrature_points  = n_quadrature_points_;
          second_n_quadrature_points = n_quadrature_points_;
 
          first_quad_map.resize(n_quadrature_points_);
          second_quad_map.resize(n_quadrature_points_);
 
          unsigned int idx = 0;
          for (const unsigned int &i : fe_values_1.quadrature_point_indices())
            for (const unsigned int &j : fe_values_2.quadrature_point_indices())
              {
                first_quad_map[idx]  = i;
                second_quad_map[idx] = j;
                ++idx;
              }
          break;
        }
      case QuadratureCouplingType::unrolled:
        {
          Assert(fe_values_1.get_quadrature_points().size() ==
                   fe_values_2.get_quadrature_points().size(),
                 ExcMessage("The two FEValuesBase objects must have the same "
                            "number of quadrature points"));
 
          n_quadrature_points_ = fe_values_1.get_quadrature_points().size();
 
          first_n_quadrature_points  = n_quadrature_points_;
          second_n_quadrature_points = n_quadrature_points_;
 
          first_quad_map.clear();
          second_quad_map.clear();
 
          break;
        }
      case QuadratureCouplingType::matching:
        {
          const auto &quadrature_points_1 = fe_values_1.get_quadrature_points();
          const auto &quadrature_points_2 = fe_values_2.get_quadrature_points();
 
          Assert(quadrature_points_1.size() == quadrature_points_2.size(),
                 ExcMessage("The two FEValuesBase objects must have the same "
                            "number of quadrature points"));
 
          for (const unsigned int &i : fe_values_1.quadrature_point_indices())
            {
              Assert(quadrature_points_1[i].distance(quadrature_points_2[i]) <
                       tol,
                     ExcMessage(
                       "The two FEValuesBase objects must have the same "
                       "quadrature points"));
            }
 
          n_quadrature_points_ = fe_values_1.get_quadrature_points().size();
 
          first_n_quadrature_points  = n_quadrature_points_;
          second_n_quadrature_points = n_quadrature_points_;
 
          first_quad_map.clear();
          second_quad_map.clear();
 
          break;
        }
      case QuadratureCouplingType::reorder:
        {
          const auto &quadrature_points_1 = fe_values_1.get_quadrature_points();
          const auto &quadrature_points_2 = fe_values_2.get_quadrature_points();
 
          Assert(quadrature_points_1.size() == quadrature_points_2.size(),
                 ExcMessage("The two FEValuesBase objects must have the same "
                            "number of quadrature points"));
 
          n_quadrature_points_ = fe_values_1.get_quadrature_points().size();
 
          first_n_quadrature_points  = n_quadrature_points_;
          second_n_quadrature_points = n_quadrature_points_;
 
          // The first is the id. The second is renumbered.
          first_quad_map.clear();
          second_quad_map.resize(fe_values_1.get_quadrature_points().size());
 
          // [TODO]: Avoid quadratic complexity here
          for (const unsigned int &i : fe_values_1.quadrature_point_indices())
            {
              auto id = numbers::invalid_unsigned_int;
              for (const unsigned int &j :
                   fe_values_2.quadrature_point_indices())
                if (quadrature_points_1[i].distance(quadrature_points_2[j]) <
                    tol)
                  {
                    id                 = i;
                    second_quad_map[i] = j;
                    break;
                  }
              Assert(id != numbers::invalid_unsigned_int,
                     ExcMessage(
                       "The two FEValuesBase objects must have the same "
                       "quadrature points, even if not in the same order."));
            }
          break;
        }
      case QuadratureCouplingType::overlapping:
        {
          const auto &quadrature_points_1 = fe_values_1.get_quadrature_points();
          const auto &quadrature_points_2 = fe_values_2.get_quadrature_points();
 
          // [TODO]: Avoid quadratic complexity here
          for (const unsigned int &i : fe_values_1.quadrature_point_indices())
            {
              for (const unsigned int &j :
                   fe_values_2.quadrature_point_indices())
                if (quadrature_points_1[i].distance(quadrature_points_2[j]) <
                    tol)
                  {
                    first_quad_map.emplace_back(i);
                    second_quad_map.emplace_back(j);
                    break;
                  }
            }
          n_quadrature_points_ = first_quad_map.size();
 
          first_n_quadrature_points  = n_quadrature_points_;
          second_n_quadrature_points = n_quadrature_points_;
          break;
        }
      default:
        Assert(false, ExcInternalError());
    }
 
  first_renumbering_data = std::make_unique<FEValuesViews::RenumberingData>(
    first_n_inner_dofs,
    first_n_inner_quadrature_points,
    first_dofs_map,
    first_quad_map);
 
  second_renumbering_data = std::make_unique<FEValuesViews::RenumberingData>(
    second_n_inner_dofs,
    second_n_inner_quadrature_points,
    second_dofs_map,
    second_quad_map);
}
 
 
template <int dim1, int dim2, int spacedim>
inline double
FECouplingValues<dim1, dim2, spacedim>::JxW(const unsigned int q) const
{
  AssertIndexRange(q, n_quadrature_points_);
 
  const auto first_q = first_renumbering_data->quadrature_renumbering.empty() ?
                         q :
                         first_renumbering_data->quadrature_renumbering[q];
 
  if (quadrature_coupling_type == QuadratureCouplingType::tensor_product ||
      quadrature_coupling_type == QuadratureCouplingType::unrolled)
    {
      const auto second_q =
        second_renumbering_data->quadrature_renumbering.empty() ?
          q :
          second_renumbering_data->quadrature_renumbering[q];
      return first_fe_values->JxW(first_q) * second_fe_values->JxW(second_q);
    }
  else
    return first_fe_values->JxW(first_q);
}
 
 
 
template <int dim1, int dim2, int spacedim>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FECouplingValues<dim1, dim2, spacedim>::quadrature_point_indices() const
{
  return {0U, n_quadrature_points_};
}
 
 
 
template <int dim1, int dim2, int spacedim>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FECouplingValues<dim1, dim2, spacedim>::coupling_dof_indices() const
{
  AssertThrow(n_coupling_dofs_ != numbers::invalid_unsigned_int,
              ExcMessage(
                "Dofs are independent. You cannot ask for coupling dofs."));
  return {0U, n_coupling_dofs_};
}
 
 
 
template <int dim1, int dim2, int spacedim>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FECouplingValues<dim1, dim2, spacedim>::first_dof_indices() const
{
  return {0U, n_first_dofs()};
}
 
 
 
template <int dim1, int dim2, int spacedim>
inline std_cxx20::ranges::iota_view<unsigned int, unsigned int>
FECouplingValues<dim1, dim2, spacedim>::second_dof_indices() const
{
  return {0U, n_second_dofs()};
}
 
 
 
template <int dim1, int dim2, int spacedim>
inline unsigned int
FECouplingValues<dim1, dim2, spacedim>::n_coupling_dofs() const
{
  AssertThrow(n_coupling_dofs_ != numbers::invalid_unsigned_int,
              ExcMessage(
                "Dofs are independent. You cannot ask for coupling dofs."));
  return n_coupling_dofs_;
}
 
 
 
template <int dim1, int dim2, int spacedim>
inline unsigned int
FECouplingValues<dim1, dim2, spacedim>::n_first_dofs() const
{
  return first_renumbering_data->n_dofs;
}
 
 
 
template <int dim1, int dim2, int spacedim>
inline unsigned int
FECouplingValues<dim1, dim2, spacedim>::n_second_dofs() const
{
  return second_renumbering_data->n_dofs;
}
 
 
 
template <int dim1, int dim2, int spacedim>
inline unsigned int
FECouplingValues<dim1, dim2, spacedim>::n_quadrature_points() const
{
  return n_quadrature_points_;
}
 
 
 
template <int dim1, int dim2, int spacedim>
std::pair<Point<spacedim>, Point<spacedim>>
FECouplingValues<dim1, dim2, spacedim>::quadrature_point(
  const unsigned int quadrature_point) const
{
  AssertIndexRange(quadrature_point, n_quadrature_points_);
  const auto first_q = first_fe_values->quadrature_point(
    first_renumbering_data->quadrature_renumbering.empty() ?
      quadrature_point :
      first_renumbering_data->quadrature_renumbering[quadrature_point]);
 
  const auto second_q = second_fe_values->quadrature_point(
    second_renumbering_data->quadrature_renumbering.empty() ?
      quadrature_point :
      second_renumbering_data->quadrature_renumbering[quadrature_point]);
 
  return {first_q, second_q};
}
 
 
 
template <int dim1, int dim2, int spacedim>
std::pair<unsigned int, unsigned int>
FECouplingValues<dim1, dim2, spacedim>::
  coupling_quadrature_to_quadrature_indices(
    const unsigned int quadrature_point) const
{
  AssertIndexRange(quadrature_point, n_quadrature_points_);
  const auto first_id =
    first_renumbering_data->quadrature_renumbering.empty() ?
      quadrature_point :
      first_renumbering_data->quadrature_renumbering[quadrature_point];
 
  const auto second_id =
    second_renumbering_data->quadrature_renumbering.empty() ?
      quadrature_point :
      second_renumbering_data->quadrature_renumbering[quadrature_point];
  return std::make_pair(first_id, second_id);
}
 
 
 
template <int dim1, int dim2, int spacedim>
std::pair<unsigned int, unsigned int>
FECouplingValues<dim1, dim2, spacedim>::coupling_dof_to_dof_indices(
  const unsigned int coupling_dof_index) const
{
  AssertIndexRange(coupling_dof_index, n_coupling_dofs_);
  const auto first_id =
    first_renumbering_data->dof_renumbering.empty() ?
      coupling_dof_index :
      first_renumbering_data->dof_renumbering[coupling_dof_index];
 
  const auto second_id =
    second_renumbering_data->dof_renumbering.empty() ?
      coupling_dof_index :
      second_renumbering_data->dof_renumbering[coupling_dof_index];
  return std::make_pair(first_id, second_id);
}
 
 
 
template <int dim1, int dim2, int spacedim>
template <typename Extractor>
inline FEValuesExtractors::FirstCoupling<Extractor>
FECouplingValues<dim1, dim2, spacedim>::get_first_extractor(
  const Extractor &extractor) const
{
  return FEValuesExtractors::FirstCoupling<Extractor>(extractor);
}
 
 
 
template <int dim1, int dim2, int spacedim>
template <typename Extractor>
inline FEValuesExtractors::SecondCoupling<Extractor>
FECouplingValues<dim1, dim2, spacedim>::get_second_extractor(
  const Extractor &extractor) const
{
  return FEValuesExtractors::SecondCoupling<Extractor>(extractor);
}
 
 
 
template <int dim1, int dim2, int spacedim>
template <typename Extractor>
inline const FEValuesViews::RenumberedView<
  FEValuesViews::View<dim1, spacedim, Extractor>>
FECouplingValues<dim1, dim2, spacedim>::operator[](
  const FEValuesExtractors::FirstCoupling<Extractor> &extractor) const
{
  return FEValuesViews::RenumberedView<
    typename FEValuesViews::View<dim1, spacedim, Extractor>>(
    (*first_fe_values)[extractor.extractor], *first_renumbering_data);
}
 
 
 
template <int dim1, int dim2, int spacedim>
template <typename Extractor>
inline const FEValuesViews::RenumberedView<
  typename FEValuesViews::View<dim2, spacedim, Extractor>>
FECouplingValues<dim1, dim2, spacedim>::operator[](
  const FEValuesExtractors::SecondCoupling<Extractor> &extractor) const
{
  return FEValuesViews::RenumberedView<
    typename FEValuesViews::View<dim2, spacedim, Extractor>>(
    (*second_fe_values)[extractor.extractor], *second_renumbering_data);
}
 
#endif // DOXYGEN
 
DEAL_II_NAMESPACE_CLOSE
 
#endif