]> https://gitweb.dealii.org/ - dealii.git/commitdiff
FECouplingValues. 15773/head
authorLuca Heltai <luca.heltai@sissa.it>
Thu, 20 Jul 2023 22:54:48 +0000 (00:54 +0200)
committerLuca Heltai <luca.heltai@sissa.it>
Thu, 11 Apr 2024 13:20:44 +0000 (15:20 +0200)
14 files changed:
doc/news/changes/minor/20230721LucaHeltai [new file with mode: 0644]
include/deal.II/fe/fe_coupling_values.h
include/deal.II/fe/fe_values_extractors.h
tests/fecoupling/CMakeLists.txt [new file with mode: 0644]
tests/fecoupling/fe_coupling_01.cc [new file with mode: 0644]
tests/fecoupling/fe_coupling_01.output [new file with mode: 0644]
tests/fecoupling/fe_coupling_02.cc [new file with mode: 0644]
tests/fecoupling/fe_coupling_02.output [new file with mode: 0644]
tests/fecoupling/fe_coupling_03.cc [new file with mode: 0644]
tests/fecoupling/fe_coupling_03.output [new file with mode: 0644]
tests/fecoupling/fe_coupling_04.cc [new file with mode: 0644]
tests/fecoupling/fe_coupling_04.output [new file with mode: 0644]
tests/fecoupling/fe_coupling_05.cc [new file with mode: 0644]
tests/fecoupling/fe_coupling_05.output [new file with mode: 0644]

diff --git a/doc/news/changes/minor/20230721LucaHeltai b/doc/news/changes/minor/20230721LucaHeltai
new file mode 100644 (file)
index 0000000..9caf84b
--- /dev/null
@@ -0,0 +1,4 @@
+New: Created a new class FECouplingValues, that helps with assembling of coupled 
+finite element methods across different dimensions or different grids.
+<br>
+(Luca Heltai, 2023/07/21)
index c6ab804fecdf0d0ab7be2482d72c4a7b1cf2f1bd..2ed108beec311dddba5a41dc9e9b3cc210593a93 100644 (file)
@@ -20,6 +20,7 @@
 
 #include <deal.II/algorithms/general_data_storage.h>
 
+#include <deal.II/base/std_cxx20/iota_view.h>
 #include <deal.II/base/subscriptor.h>
 #include <deal.II/base/thread_local_storage.h>
 #include <deal.II/base/utilities.h>
@@ -406,6 +407,664 @@ namespace FEValuesViews
 } // namespace FEValuesViews
 
 
+/**
+ * Quadrature coupling options when assembling quadrature formulas for double
+ * integrals.
+ *
+ * When computing the approximation of double integrals of the form
+ *
+ * \f[
+ * \int_{T_1} \int{T_2} K(x_1, x_2) f(x_1) g(x_2) dT_1 dT_2,
+ * \f]
+ *
+ * where $T_1$ and $T_2$ are two arbitrary sets (cells, faces, edges, or any
+ * combination thereof), and $K$ is a (possibly singular) coupling kernel, one
+ * needs to combine quadrature formulas from two different FEValuesBase objects.
+ *
+ * This enum class provides a way to specify how the quadrature points and
+ * weights should be combined. In general, the two FEValuesBase objects provide
+ * different quadrature rules, and these can be interpreted in different ways,
+ * depending on the kernel function that is being integrated, and on how the two
+ * quadrature rules were constructed.
+ *
+ * This enum is used in the constructor of FECouplingValues to specify how to
+ * interpret and manipulate the quadrature points and weights of the two
+ * FEValuesBase objects.
+ */
+enum class QuadratureCouplingType
+{
+  /**
+   * The FEValuesBase objects provide different quadrature rules, and the
+   * resulting quadrature points and weights should be constructed as the tensor
+   * product of the two quadrature rules. This is generally used for non-local
+   * and non-singular kernels, where the quadrature points of the two
+   * FEValuesBase objects are independent of each other.
+   */
+  tensor_product,
+
+  /**
+   * Both FEValuesBase objects provide the same number of (generally different)
+   * quadrature points, that should be used as is to implement the double
+   * integration. This is equivalent to rewriting double integrals as a single
+   * sum over unrolled indices, and it is useful, for example, when performing
+   * integration of singular kernels, which require special quadrature rules
+   * both on the inside and on the outside integral. An exception is thrown if
+   * the two FEValuesBase objects do not provide the same number of quadrature
+   * points, but otherwise the two sets of points and weights can be arbitrary.
+   */
+  unrolled,
+
+  /**
+   * The FEValuesBase objects provide the same quadrature rule over the same
+   * set, with the same orientation. This is similar to the unrolled case, in
+   * the sense that the quadrature formulas are used as is, but both objects are
+   * assumed to return the same quadrature points and weights. This is useful
+   * when integrating over faces of neighboring cells, and when the reordering
+   * of the quadrature points is known to match a priori. An exception is thrown
+   * if the quadrature points and weights of the two FEValuesBase objects do not
+   * match one-to-one.
+   */
+  matching,
+
+  /**
+   * The FEValuesBase objects provide the same quadrature rule over the same
+   * set, but with possibly different orientations. The quadrature points are
+   * reordered internally, so that the resulting quadrature points and weights
+   * are the same as in the matching case. An exception is thrown if the
+   * quadrature points and weights of the two FEValuesBase cannot be reordered
+   * to match one-to-one.
+   */
+  reorder,
+
+  /**
+   * The FEValuesBase objects provide partially overlapping quadrature rules
+   * over two intersecting sets, with possibly different orientations. The
+   * quadrature points are reordered and matched internally, so that the
+   * resulting quadrature points and weights are only the ones that match
+   * between the two objects. If no overlapping points and weights can be found,
+   * an empty set is used to integrate, and the resulting integral is zero. No
+   * exception is thrown in this case.
+   */
+  overlapping,
+};
+
+/**
+ * DoF coupling options when assembling double integrals.
+ *
+ * When computing the approximation of double integrals of the form
+ *
+ * \f[
+ * \int_{T_1} \int{T_2} K(x_1, x_2) v_i(x_1) w_j(x_2) dT_1 dT_2,
+ * \f]
+ *
+ * where $T_1$ and $T_2$ are two arbitrary sets (cells, faces, edges, or any
+ * combination thereof), and $K$ is a (possibly singular) coupling kernel, one
+ * may want to combine degrees from two different FEValuesBase objects (i.e.,
+ * basis functions $v_i$ and $w_j$ in the examples above).
+ *
+ * This enum class provides a way to specify how the degrees of freedom should
+ * be combined. There are two cases of interest:
+ *
+ * 1. the two FEValuesBase objects refer to different DoFHandlers
+ * 2. the two FEValuesBase objects refer to the same DoFHandler
+ *
+ * In the first case, one usually treats the two sets of degrees of freedom as
+ * independent of each other, and the resulting matrix is generally rectangular.
+ *
+ * In the second case, one may choose to treat the two sets of degrees of
+ * freedom either as independent or to group them together. A similar approach
+ * is used in the FEInterfaceValues class, where the degrees of freedom of the
+ * two FEValuesBase objects are grouped together, in a contiguous way, so that
+ * the resulting basis functions are interpreted in the following way:
+ *
+ * \f[
+ * \phi_{1,i}(x) = \begin{cases} v_i(x) & \text{ if } i \in [0,n_l) \\
+ * 0 & \text{ if ) i \in [n_1, n_1+n_2] \end{cases},\quad \phi_{1,i}(x) =
+ * \begin{cases} 0(x) & \text{ if } i \in [0,n_1) \\
+ * w_{i-n_1}(x) & \text{ if ) i \in [n_1, n_1+n_2] \end{cases},
+ * \f]
+ *
+ * where $phi_{1,i}$ is the first basis function with index $i$ and $n_{1,2}$
+ * are the number of local dofs on the first and second FEValuesBase objects.
+ *
+ * This enum is used in the constructor of FECouplingValues to specify how to
+ * interpret and manipulate the local dof indices of the two FEValuesBase
+ * objects.
+ */
+enum class DoFCouplingType
+{
+  /**
+   * The FEValuesBase objects may have different dof indices, possibly indexing
+   * different DoFHandler objects, and we are interested in assembling a
+   * generally rectangular matrix, where there is no relationship between the
+   * two index spaces.
+   */
+  independent,
+
+  /**
+   * The FEValuesBase objects may have different dof indices, but they both
+   * index the same DoFHandler objects, and we are interested in assembling them
+   * all together in a single matrix. In this coupling type, the DoF indices are
+   * grouped together, one after the other, first the first dof indices, and
+   * then the second dof indices. This is useful when one wants to assemble four
+   * matrices at the same time, corresponding to the four possible combinations
+   * of the two FEValuesBase objects, (i.e., first coupled with first, first
+   * coupled with second, second coupled with first, and second coupled with
+   * second) and then sum them together to obtain the final system. This is
+   * similar to what is done in the FEInterfaceValues class.
+   */
+  contiguous,
+};
+
+/**
+ * FECouplingValues is a class that facilitates the integration of finite
+ * element data between two different finite element objects, possibly living on
+ * different grids, and with possibly different topological dimensions (i.e.,
+ * cells, faces, edges, and any combination thereof).
+ *
+ * This class provides a way to simplify the implementation of the following
+ * abstract operation:
+ *
+ * \f[
+ * \int_{T_1} \int{T_2} K(x_1, x_2) \phi^1_i(x_1) \phi^2_j(x_2) dT_1 dT_2
+ * \f]
+ *
+ * for three different types of Kernels $K$:
+ * - $K(x_1, x_2)$ is a non-singular Kernel function, for example, it is a
+ *   function of positive powers $\alpha$ of the distance between the quadrature
+ *   points $|x_1-x_2|^\alpha$;
+ * - $K(x_1, x_2)$ is a singular Kernel function, for example, it is a function
+ *   of negative powers $\alpha$ of the distance between the quadrature points
+ *   $|x_1-x_2|^\alpha$;
+ * - $K(x_1, x_2)$ is a Dirac delta distribution $\delta(x_1-x_2)$, such that
+ *   the integral above is actually a single integral over the intersection of
+ *   the two sets $T_1$ and $T_2$.
+ *
+ * For the first case, one may think that the only natural way to proceed is to
+ * compute the double integral by simply nesting two loops:
+ * \f[
+ * \int_{T_1} \int{T_2} K(x_1, x_2) \phi^1_i(x_1) \phi^2_j(x_2) dT_1 dT_2
+ * \approx \sum_{q_1} \sum_{q_2} K(x_1^{q_1}, x_2^{q_2}) \phi^1_i(x_1^{q_1})
+ * \phi^2_j(x_2^{q_2}) w_1^{q_1} w_2^{q_2},
+ * \f]
+ *
+ * where $x_1^{q_1}$ and $x_2^{q_2}$ are the quadrature points in $T_1$ and
+ * $T_2$ respectively, and $w_1^{q_1}$ and $w_2^{q_2}$ are the corresponding
+ * quadrature weights.
+ *
+ * This, however is not the only way to proceed. In fact, such an integral can
+ * be rewritten as a single loop over corresponding elements of two arrays of
+ * points with the same length that can be thought of as a single quadrature
+ * rule on the set $T_1\times T_2$. For singular kernels, for example, this is
+ * often the only way to proceed, since the quadrature formula on $T_1\times
+ * T_2$ is usually not written as a tensor product quadrature formula, and one
+ * needs to build a custom quadrature formula for this purpose.
+ *
+ * This class  allows one to treat the three cases above in the same way, and to
+ * approximate the integral as follows:
+ *
+ * \f[
+ * \int_{T_1} \int{T_2} K(x_1, x_2) \phi^1_i(x_1) \phi^2_j(x_2) dT_1 dT_2
+ * \approx \sum_{i=1}^{N_q} K(x_1^{i}, x_2^{i}) \phi^1_i(x_1^{i})
+ * \phi^2_j(x_2^{i}) w_1^{i} w_2^i,
+ * \f]
+ *
+ * Since the triple of objects $(\{q\}, \{w\}, \{\phi\})$ is usually provided by
+ * a class derived from the FEValuesBase class, this is the type that the class
+ * needs at construction time. $T_1$ and $T_2$ can be two arbitrary cells,
+ * faces, or edges belonging to possibly different meshes (or to meshes with
+ * different topological dimensions), $\phi^1_i$ and $\phi^2_j$ are basis
+ * functions defined on $T_1$ and $T_2$, respectively.
+ *
+ * The case of the Dirac distribution is when $T_1$ and $T_2$
+ * correspond to the common face of two neighboring cells. In this case, this
+ * class provides a functionality which is similar to the FEInterfaceValues
+ * class, and gives you a way to access values of basis functions on the
+ * neighboring cells, as well as their gradients and Hessians, in a unified
+ * fashion, on the face.
+ *
+ * Similarly, this class can be used to couple bulk and surface meshes across
+ * the faces of the bulk mesh. In this case, the two FEValuesBase objects will
+ * have different topological dimension (i.e., one will be a cell in a
+ * co-dimension one triangulation, and the other a face of a bulk grid with
+ * co-dimension zero), and the QuadratureCouplingType argument is usually chosen
+ * to be QuadratureCouplingType::reorder, since the quadrature points of the two
+ * different FEValuesBase objects are not necessarily generated with the same
+ * ordering.
+ *
+ * The type of integral to compute is controlled by the QuadratureCouplingType
+ * argument (see the documentation of that enum class for more details), while
+ * the type degrees of freedom coupling is controlled by the DoFCouplingType
+ * argument (see the documentation of that enum class for more details).
+ *
+ * As an example usage of this class, consider the a bilinear form of the form:
+ * \f[
+ * \int_{T_1} \int{T_2} K_1(x_1, x_2) v_i(x_1) u_j(x_2) dT_1 dT_2 +
+ * \int_{T_1} \int{T_2} K_2(x_1, x_2) p_i(x_1) q_j(x_2) dT_1 dT_2
+ * \f]
+ * where the finite dimensional space has two scalar components. We indicate
+ * with $v_i$ and $p_i$ the trial functions, and with $u_j$ and
+ * $q_j$ the corresponding test functions. $K_1$ and $K_2$ are coupling kernels:
+ * such a formulation is used, for example, to write the bilinear forms of
+ * Galerkin boundary element methods.
+ *
+ * The corresponding implementation would look like the following:
+ *
+ * @code
+ * ... // double loop over cells that yields cell_1 and cell_2
+ *
+ * fe_values_1.reinit(cell_1);
+ * fe_values_2.reinit(cell_2);
+ *
+ * CouplingFEValues<dim, dim, spacedim> cfv(fe_values1, fe_values2,
+ *                                          DoFCouplingType::independent,
+ *                                          QuadratureCouplingType::tensor_product);
+ *
+ * FullMatrix<double> local_matrix(fe_values1.dofs_per_cell,
+ *                                 fe_values2.dofs_per_cell);
+ *
+ * // Extractor on first cell
+ * const auto v1 = cfv.get_first_extractor(FEValuesExtractor::Scalar(0));
+ * const auto p1 = cfv.get_first_extractor(FEValuesExtractor::Scalar(1));
+ *
+ * // Extractor on second cell
+ * const auto u2 = cfv.get_second_extractor(FEValuesExtractor::Scalar(0));
+ * const auto q2 = cfv.get_second_extractor(FEValuesExtractor::Scalar(1));
+ *
+ * ...
+ *
+ * for (const unsigned int q : cfv.quadrature_point_indices()) {
+ *   const auto &[x_q,y_q] = cfv.quadrature_point(q);
+ *
+ *   for (const unsigned int i : cfv.first_dof_indices())
+ *     {
+ *       const auto &v_i = cfv[v1].value(i, q);
+ *       const auto &p_i = cfv[p1].value(i, q);
+ *
+ *       for (const unsigned int j : cfv.second_dof_indices())
+ *         {
+ *           const auto &u_j = cfv[u2].value(j, q);
+ *           const auto &q_j = cfv[q2].value(j, q);
+ *
+ *           local_matrix(i, j) += (K1(x_q, y_q) * v_i * u_j +
+ *                                  K2(x_q, y_q) * p_i * q_j) *
+ *                                 cfv[0].JxW(q) *
+ *                                 cfv[1].JxW(q);
+ *         }
+ *     }
+ * }
+ * @endcode
+ *
+ * In the above loop, the quadrature points of the two FEValuesBase objects are
+ * grouped together in a single loop, while the dof indices are kept separate.
+ *
+ * Internally, this class provides an abstraction to organize coupling terms
+ * between two arbitrary FEValuesBase objects, and provides a unified way to
+ * access the values of the two basis, and of the two sets of quadrature points
+ * and weights. The underlying data is stored in the two FEValuesBase objects,
+ * and this class provides access to two RenumberedView objects, which can be
+ * obtained by calling the `operator[]` method of this class with the arguments
+ * constructed by calling the get_first_extractor() and get_second_extractor()
+ * methods of this class with a standard FEValuesExtractors object, i.e.,
+ *
+ * @code
+ * // extract the first scalar component of the basis
+ * FEValuesExtractor scalar(0);
+ * ...
+ *
+ * FECouplingValues<dim> fe_coupling(fev1, fev1);
+ *
+ * // Extractors for the two FEValuesBase objects are returned by the
+ * // get_first_extractor() and
+ * // get_second_extractor() methods
+ *
+ * const auto first = fe_coupling.get_first_extractor(scalar);
+ * const auto second = fe_coupling.get_second_extractor(scalar);
+ *
+ * // Now we can use the augmented extractors to access the values of the two
+ * // FEValuesBase objects
+ * const auto & first_vi = fe_coupling[first].value(i, q);
+ * const auto & second_vi = fe_coupling[second].value(i, q);
+ * @endcode
+ */
+template <int dim1, int dim2 = dim1, int spacedim = dim1>
+class FECouplingValues
+{
+public:
+  /**
+   * Construct the FECouplingValues in an invalid state. Before you can use this
+   * object, you must call the reinit() function.
+   */
+  FECouplingValues();
+
+  /**
+   * Construct the FECouplingValues with two arbitrary FEValuesBase objects.
+   * This class assumes that the FEValuesBase objects that are given at
+   * construction time are initialized and ready to use (i.e., that you have
+   * called the reinit() function on them before calling this constructor).
+   *
+   * Notice that the actual renumbering of the degrees of freedom and quadrature
+   * points is done at construction time, or upon calling the reinit() function.
+   * If you change the underlying FEValuesBase objects after construction, you
+   * must call the reinit() function to update the renumbering. This may or may
+   * not be necessary, depending on the type of coupling that you are using.
+   *
+   * This really depends on the application and on the specific type of
+   * coupling. For example, for volume/volume coupling, i.e., for operators with
+   * non-local and non-singular kernels of type
+   * \f[
+   * \int_K \int_T f(\phi_i(x)-\phi_j(y), x-y) dx dy
+   * \f]
+   * you may initialize FECouplingValues once, and just reinit the underlying
+   * FEValuesBase objects on different cells K and T, without the need to
+   * recompute the coupling (i.e., the numbering is always the same, and nothing
+   * differs from what happened in the first call).
+   *
+   * For cell/surface coupling, the same cell may couple with different faces,
+   * so the renumbering must be really computed from scratch for each pair of
+   * FEValuesBase objects, so reinitializing the underlying cells and faces will
+   * make the renumbering itself invalid, and FECouplingValues must be
+   * reinitialized (o constructed from scratch) after calling
+   * `fe_values_1.reinit(K)` and `fe_values_1.reinit(T)`.
+   *
+   * @param fe_values_1 The first FEValuesBase object.
+   * @param fe_values_2 The second FEValuesBase object.
+   * @param dof_coupling_type The type of dof coupling to use.
+   * @param quadrature_coupling_type The type of quadrature to use for the coupling.
+   */
+  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);
+
+  /**
+   * Reinitialize the FECouplingValues with two arbitrary FEValuesBase objects.
+   * The FEValuesBase objects must be initialized and ready to use, i.e., you
+   * must have called the reinit() function on them before calling this method.
+   *
+   * This method computes the actual renumbering of the degrees of freedom and
+   * quadrature points. If you change the underlying FEValuesBase objects after
+   * calling this method, you mey need to call the reinit() function to update
+   * the renumbering. This may or may not be necessary, depending on the type of
+   * coupling that you are using.
+   *
+   * This really depends on the application and on the specific type of
+   * coupling. For example, for volume/volume coupling, i.e., for operators with
+   * non-local and non-singular kernels of type
+   * \f[
+   * \int_K \int_T f(\phi_i(x)-\phi_j(y), x-y) dx dy
+   * \f]
+   * you may initialize FECouplingValues once, and just reinit the underlying
+   * FEValuesBase objects on different cells K and T, without the need to
+   * recompute the coupling (i.e., the numbering is always the same, and nothing
+   * differs from what happened in the first call).
+   *
+   * For cell/surface coupling, the same cell may couple with different faces,
+   * so the renumbering must be really computed from scratch for each pair of
+   * FEValuesBase objects, so reinitializing the underlying cells and faces will
+   * make the renumbering itself invalid, and FECouplingValues must be
+   * reinitialized (o constructed from scratch) after calling
+   * `fe_values_1.reinit(K)` and `fe_values_1.reinit(T)`.
+   *
+   * @param fe_values_1 The first FEValuesBase object.
+   * @param fe_values_2 The second FEValuesBase object.
+   * @param dof_coupling_type The type of dof coupling to use.
+   * @param quadrature_coupling_type The type of quadrature to use for the coupling.
+   */
+  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);
+
+
+  /**
+   * Return a FirstCoupling object that can be used to extract values from the
+   * first FEValuesBase object.
+   */
+  template <typename Extractor>
+  FEValuesExtractors::FirstCoupling<Extractor>
+  get_first_extractor(const Extractor &extractor) const;
+
+  /**
+   * Return a SecondCoupling object that can be used to extract values from the
+   * second FEValuesBase object.
+   */
+  template <typename Extractor>
+  FEValuesExtractors::SecondCoupling<Extractor>
+  get_second_extractor(const Extractor &extractor) const;
+
+  /**
+   * Return the value of JxW at the given quadrature point.
+   *
+   * @dealiiRequiresUpdateFlags{update_JxW_values}
+   */
+  double
+  JxW(const unsigned int quadrature_point) const;
+
+  /**
+   * Return the two quadrature points in real space at the given quadrature
+   * point index, corresponding to a quadrature point in the set $T_1\times
+   * T_2$.
+   *
+   * @dealiiRequiresUpdateFlags{update_quadrature_points}
+   */
+  std::pair<Point<spacedim>, Point<spacedim>>
+  quadrature_point(const unsigned int quadrature_point) const;
+
+  /**
+   * Return an object that can be thought of as an array containing all
+   * indices from zero to `n_quadrature_points`. This allows to write code
+   * using range-based `for` loops.
+   *
+   * @see CPP11
+   */
+  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+  quadrature_point_indices() const;
+
+  /**
+   * Return an object that can be thought of as an array containing all
+   * indices from zero (inclusive) to `n_first_dofs()` (exclusive).
+   * This allows one to write code using range-based `for` loops.
+   */
+  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+  first_dof_indices() const;
+
+  /**
+   * Return an object that can be thought of as an array containing all
+   * indices from zero (inclusive) to `n_second_dofs()` (exclusive).
+   * This allows one to write code using range-based `for` loops.
+   */
+  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+  second_dof_indices() const;
+
+  /**
+   * Return an object that can be thought of as an array containing all indices
+   * from zero (inclusive) to `n_coupling_dof_indices()` (exclusive). This
+   * allows one to write code using range-based `for` loops.
+   *
+   * This function only makes sense when the dof coupling type is
+   * DoFCouplingType::contiguous, and, in this case, n_coupling_dof_indices() is
+   * the  sum of n_first_dofs() and n_second_dofs(). If
+   * DoFCouplingType::independent is used, you should call first_dof_indices()
+   * or second_dof_indices() instead, and calling this function will throw an
+   * exception.
+   */
+  std_cxx20::ranges::iota_view<unsigned int, unsigned int>
+  coupling_dof_indices() const;
+
+  /**
+   * Return the set of joint DoF indices, possibly offset by the given values.
+   *
+   * For coupled problems, the DoF indices of the two FEValuesBase objects are
+   * generally constructed from different DoFHandler objects. According to how
+   * the user wants to construct the coupling matrices, the two dof indices may
+   * map local dof indices to global dof indices starting from zero (i.e., we
+   * are assembling a single rectangular matrix, and the first dof indices run
+   * on rows of the matrix, while the second dof indices run on columns of the
+   * matrix), or the user may be wanting to map the coupling matrix to a block
+   * of a larger block matrix. In this second case the dof indices of either or
+   * both blocks may need to be shifted by a fixed amount, to fall within the
+   * right range of the target block matrix. This is usually required if the
+   * user wants to keep their solution vectors in a single block vector, where
+   * each block corresponds to a different DoFHandler object. In this case, all
+   * dofs referring to the first FEValuesBase object are offset by the value of
+   * dofs_offset_1, while all dofs referring to the second FEValuesBase object
+   * are offset by the value of dofs_offset_2, i.e., the following relationship
+   * is satisfied:
+   * @code
+   * fe_values_1.get_dof_indices(dof_indices_1);
+   * fe_values_2.get_dof_indices(dof_indices_2);
+   * auto coupling_dof_indices = fe_coupling_values.
+   *      get_coupling_dof_indices(dof_indices_1, dof_indices_2,
+   *                               dofs_offset_1, dofs_offset_2);
+   * for(const auto &i : fe_coupling_values.coupling_dof_indices()) {
+   *    auto [id1, id2] = fe_coupling_values.coupling_dof_to_dof_indices(i);
+   *    if (id1 != numbers::invalid_unsigned_int)
+   *      // always true
+   *      coupling_dof_indices[i] == dof_indices_1[id1] + dofs_offset_1;
+   *    if (id2 != numbers::invalid_unsigned_int)
+   *      // always true
+   *      coupling_dof_indices[i] == dof_indices_2[id2] + dofs_offset_2;
+   * }
+   * @endcode
+   *
+   * This method only makes sense when the dof coupling type is
+   * DoFCouplingType::contiguous. If DoFCouplingType::independent is used,
+   * calling this function will result in an exception.
+   */
+  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;
+
+  /**
+   * Convert a coupling dof index into the corresponding local DoF indices of
+   * the two FEValuesObjects. If a DoF is only active on one of the
+   * FEValuesObjects, the other index will be numbers::invalid_unsigned_int.
+   */
+  std::pair<unsigned int, unsigned int>
+  coupling_dof_to_dof_indices(const unsigned int coupling_dof_index) const;
+
+  /**
+   * Convert a quadrature index into the corresponding local quadrature indices
+   * of the two FEValuesObjects. Both indices are guaranteed to be valid within
+   * the corresponding FEValuesObject.
+   */
+  std::pair<unsigned int, unsigned int>
+  coupling_quadrature_to_quadrature_indices(
+    const unsigned int quadrature_point) const;
+
+  /**
+   * @name Extractors Methods to extract individual components
+   * @{
+   */
+
+  /**
+   * Create a combined view of the first FECouplingValues object that represents
+   * a view of the possibly vector-valued finite element. The concept of views
+   * is explained in the documentation of the namespace FEValuesViews.
+   */
+  template <typename Extractor>
+  const FEValuesViews::RenumberedView<
+    FEValuesViews::View<dim1, spacedim, Extractor>>
+  operator[](
+    const FEValuesExtractors::FirstCoupling<Extractor> &extractor) const;
+
+  /**
+   * Create a combined view of the second FECouplingValues object that
+   * represents a view of the possibly vector-valued finite element. The concept
+   * of views is explained in the documentation of the namespace FEValuesViews.
+   */
+  template <typename Extractor>
+  const FEValuesViews::RenumberedView<
+    FEValuesViews::View<dim2, spacedim, Extractor>>
+  operator[](
+    const FEValuesExtractors::SecondCoupling<Extractor> &extractor) const;
+
+  /**
+   * @}
+   */
+
+  /**
+   * Return the number of coupling DoF indices. If DoFCouplingType::independent
+   * is used, this function returns numbers::invalid_unsigned_int, otherwise it
+   * returns the sum of n_first_dofs() and n_second_dofs().
+   */
+  unsigned int
+  n_coupling_dofs() const;
+
+  /**
+   * Return the number of first DoF indices. This generally coincides with
+   * n_dofs_per_cell of the first FEValuesBase object.
+   */
+  unsigned int
+  n_first_dofs() const;
+
+  /**
+   * Return the number of second DoF indices. This generally coincides with
+   * n_dofs_per_cell of the second FEValuesBase object.
+   */
+  unsigned int
+  n_second_dofs() const;
+
+  /**
+   * Return the number of quadrature points.
+   */
+  unsigned int
+  n_quadrature_points() const;
+
+private:
+  /**
+   * The dof coupling type used by this object.
+   */
+  DoFCouplingType dof_coupling_type;
+
+  /**
+   * The quadrature coupling type used by this object.
+   */
+  QuadratureCouplingType quadrature_coupling_type;
+
+  /**
+   * Pointer to first FEValuesBase object.
+   */
+  SmartPointer<const FEValuesBase<dim1, spacedim>> first_fe_values;
+
+  /**
+   * Pointer to second FEValuesBase object.
+   */
+  SmartPointer<const FEValuesBase<dim2, spacedim>> second_fe_values;
+
+  /**
+   * Renumbering data for the first FEValuesBase object.
+   */
+  std::unique_ptr<const FEValuesViews::RenumberingData> first_renumbering_data;
+
+  /**
+   * Renumbering data for the second FEValuesBase object.
+   */
+  std::unique_ptr<const FEValuesViews::RenumberingData> second_renumbering_data;
+
+  /**
+   * Number of quadrature points.
+   */
+  unsigned int n_quadrature_points_;
+
+  /**
+   * Number of coupling DoF indices. If DoFCouplingType::independent is used,
+   * this is numbers::invalid_unsigned_int, while it is n_first_dofs() +
+   * n_second_dofs() in the the case of DoFCouplingType::contiguous.
+   */
+  unsigned int n_coupling_dofs_;
+};
+
+
 #ifndef DOXYGEN
 
 
@@ -660,6 +1319,505 @@ namespace FEValuesViews
   }
 } // 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
index 2f73c16eb5dac84796e9bd16e46aa293e1e6b0f6..ef109bca01f8232d688059803880178812d2b5f3 100644 (file)
@@ -268,6 +268,42 @@ namespace FEValuesExtractors
     std::string
     get_name() const;
   };
+
+  /**
+   * Helper struct to associate an extractor to the first FEValuesBase
+   * sub-object of a FECouplingValues object.
+   */
+  template <typename Extractor>
+  struct FirstCoupling
+  {
+    /**
+     * Construct a new FirstCoupling object with the given extractor.
+     */
+    FirstCoupling(const Extractor &extractor);
+
+    /**
+     * The actual extractor object.
+     */
+    const Extractor extractor;
+  };
+
+  /**
+   * Helper struct to associate an extractor to the second FEValuesBase
+   * sub-object of a FECouplingValues object.
+   */
+  template <typename Extractor>
+  struct SecondCoupling
+  {
+    /**
+     * Construct a new SecondCoupling object with the given extractor.
+     */
+    SecondCoupling(const Extractor &extractor);
+
+    /**
+     * The actual extractor object.
+     */
+    const Extractor extractor;
+  };
 } // namespace FEValuesExtractors
 
 
@@ -320,9 +356,20 @@ namespace FEValuesExtractors
   inline Tensor<rank>::Tensor(const unsigned int first_tensor_component)
     : first_tensor_component(first_tensor_component)
   {}
-} // namespace FEValuesExtractors
 
 
+  template <typename Extractor>
+  inline FirstCoupling<Extractor>::FirstCoupling(const Extractor &extractor)
+    : extractor(extractor)
+  {}
+
+
+  template <typename Extractor>
+  inline SecondCoupling<Extractor>::SecondCoupling(const Extractor &extractor)
+    : extractor(extractor)
+  {}
+} // namespace FEValuesExtractors
+
 
 DEAL_II_NAMESPACE_CLOSE
 
diff --git a/tests/fecoupling/CMakeLists.txt b/tests/fecoupling/CMakeLists.txt
new file mode 100644 (file)
index 0000000..db9c1ad
--- /dev/null
@@ -0,0 +1,4 @@
+cmake_minimum_required(VERSION 3.13.4)
+include(../scripts/setup_testsubproject.cmake)
+project(testsuite CXX)
+deal_ii_pickup_tests()
diff --git a/tests/fecoupling/fe_coupling_01.cc b/tests/fecoupling/fe_coupling_01.cc
new file mode 100644 (file)
index 0000000..8ba61b7
--- /dev/null
@@ -0,0 +1,224 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test basic properties of FECouplingValues
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_coupling_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../test_grids.h"
+
+
+template <typename FCV>
+void
+inspect_fcv(const FCV                                  &fcv,
+            const std::vector<types::global_dof_index> &v1,
+            const std::vector<types::global_dof_index> &v2)
+{
+  deallog << "n_coupling_dofs: " << fcv.n_coupling_dofs() << std::endl;
+
+  auto indices = fcv.get_coupling_dof_indices(v1, v2);
+  Assert(indices.size() == fcv.n_coupling_dofs(), ExcInternalError());
+
+  deallog << "coupling_dof_indices: ";
+  for (const auto &i : indices)
+    deallog << i << ' ';
+  deallog << std::endl;
+
+  unsigned int idx = 0;
+  for (const auto &v : indices)
+    {
+      const auto pair = fcv.coupling_dof_to_dof_indices(idx);
+
+      deallog << "  index " << idx << " global_dof_index: " << v
+              << ", index 1: " << static_cast<int>(pair.first)
+              << ", index 2: " << static_cast<int>(pair.second) << std::endl;
+      ++idx;
+    }
+
+  deallog << "n_quadrature_points: " << fcv.n_quadrature_points() << std::endl;
+
+  for (const auto &q : fcv.quadrature_point_indices())
+    {
+      const auto &[id1, id2] = fcv.coupling_quadrature_to_quadrature_indices(q);
+      const auto &[q1, q2]   = fcv.quadrature_point(q);
+      deallog << "  q index " << q << ", index 1: " << id1
+              << ", index 2: " << id2 << ", q1: " << q1 << ", q2: " << q2
+              << std::endl;
+    }
+
+  deallog << std::endl;
+}
+
+
+
+template <int dim>
+void
+test()
+{
+  Triangulation<dim> tria;
+  TestGrids::hyper_line(tria, 2);
+
+  DoFHandler<dim> dofh(tria);
+  FE_Q<dim>       fe(1);
+  dofh.distribute_dofs(fe);
+
+  MappingQ<dim> mapping(1);
+  UpdateFlags   update_flags = update_quadrature_points | update_JxW_values;
+
+  FEValues<dim> fv1(mapping, fe, QGauss<dim>(fe.degree + 1), update_flags);
+  FEValues<dim> fv2(mapping, fe, QGauss<dim>(fe.degree + 1), update_flags);
+
+  FEValues<dim> fv3(mapping, fe, QTrapezoid<dim>(), update_flags);
+  FEValues<dim> fv4(mapping, fe, QTrapezoid<dim>(), update_flags);
+
+  FEFaceValues<dim> fv5(mapping, fe, QTrapezoid<dim - 1>(), update_flags);
+  FEFaceValues<dim> fv6(mapping, fe, QTrapezoid<dim - 1>(), update_flags);
+
+  std::vector<types::global_dof_index> v1(fe.dofs_per_cell);
+  std::vector<types::global_dof_index> v2(fe.dofs_per_cell);
+
+  auto cell1 = dofh.begin();
+  fv1.reinit(cell1);
+  cell1->get_dof_indices(v1);
+
+  unsigned int boundary_face = numbers::invalid_unsigned_int;
+  unsigned int internal_face = numbers::invalid_unsigned_int;
+
+  for (const unsigned int f : cell1->face_indices())
+    {
+      if (cell1->at_boundary(f))
+        {
+          boundary_face = f;
+        }
+      else if (!cell1->at_boundary(f))
+        {
+          internal_face = f;
+        }
+      if (boundary_face != numbers::invalid_unsigned_int &&
+          internal_face != numbers::invalid_unsigned_int)
+        break;
+    }
+
+  auto cell2 = cell1->neighbor(internal_face);
+  fv2.reinit(cell2);
+  cell2->get_dof_indices(v2);
+
+  {
+    deallog << "** coupling between cell1 and cell2 (tensor) **" << std::endl;
+    FECouplingValues<dim> fcv(fv1, fv2, DoFCouplingType::contiguous);
+    inspect_fcv(fcv, v1, v2);
+  }
+  {
+    deallog << "** coupling between cell1 and cell2 (unrolled) **" << std::endl;
+    FECouplingValues<dim> fcv(fv1,
+                              fv2,
+                              DoFCouplingType::contiguous,
+                              QuadratureCouplingType::unrolled);
+    inspect_fcv(fcv, v1, v2);
+  }
+  {
+    deallog << "** coupling between cell1 and cell2 Trapez (overlapping) **"
+            << std::endl;
+    fv3.reinit(cell1);
+    fv4.reinit(cell2);
+    FECouplingValues<dim> fcv(fv3,
+                              fv4,
+                              DoFCouplingType::contiguous,
+                              QuadratureCouplingType::overlapping);
+    inspect_fcv(fcv, v1, v2);
+  }
+  {
+    deallog
+      << "** coupling between cell1 and boundary face, trapez (overlapping) **"
+      << std::endl;
+    fv3.reinit(cell1);
+    fv5.reinit(cell1, boundary_face);
+    FECouplingValues<dim> fcv(fv3,
+                              fv5,
+                              DoFCouplingType::contiguous,
+                              QuadratureCouplingType::overlapping);
+    inspect_fcv(fcv, v1, v1);
+  }
+  {
+    deallog
+      << "** coupling between cell1 and cell1->internal_face, trapez (overlapping) **"
+      << std::endl;
+    fv3.reinit(cell1);
+    fv5.reinit(cell1, internal_face);
+    FECouplingValues<dim> fcv(fv3,
+                              fv5,
+                              DoFCouplingType::contiguous,
+                              QuadratureCouplingType::overlapping);
+    inspect_fcv(fcv, v1, v2);
+  }
+  {
+    deallog
+      << "** coupling between cell1 and cell2->internal_face, trapez (overlapping) **"
+      << std::endl;
+    fv3.reinit(cell1);
+    fv5.reinit(cell2, cell1->neighbor_of_neighbor(internal_face));
+    FECouplingValues<dim> fcv(fv3,
+                              fv5,
+                              DoFCouplingType::contiguous,
+                              QuadratureCouplingType::overlapping);
+    inspect_fcv(fcv, v1, v2);
+  }
+  {
+    deallog
+      << "** coupling between cell1 and cell2->other_face, trapez (overlapping) **"
+      << std::endl;
+    fv3.reinit(cell1);
+    fv5.reinit(cell2, internal_face);
+    FECouplingValues<dim> fcv(fv3,
+                              fv5,
+                              DoFCouplingType::contiguous,
+                              QuadratureCouplingType::overlapping);
+    inspect_fcv(fcv, v1, v2);
+  }
+  {
+    deallog
+      << "** coupling between cell1->face(0) and cell1->face(2), trapez (tensor) **"
+      << std::endl;
+    fv5.reinit(cell1, 0);
+    fv6.reinit(cell1, 2);
+    FECouplingValues<dim> fcv(fv5,
+                              fv6,
+                              DoFCouplingType::contiguous,
+                              QuadratureCouplingType::tensor_product);
+    inspect_fcv(fcv, v1, v1);
+  }
+}
+
+int
+main()
+{
+  initlog();
+  test<2>();
+}
diff --git a/tests/fecoupling/fe_coupling_01.output b/tests/fecoupling/fe_coupling_01.output
new file mode 100644 (file)
index 0000000..1e45d52
--- /dev/null
@@ -0,0 +1,137 @@
+
+DEAL::** coupling between cell1 and cell2 (tensor) **
+DEAL::n_coupling_dofs: 8
+DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 
+DEAL::  index 0 global_dof_index: 0, index 1: 0, index 2: -1
+DEAL::  index 1 global_dof_index: 1, index 1: 1, index 2: -1
+DEAL::  index 2 global_dof_index: 2, index 1: 2, index 2: -1
+DEAL::  index 3 global_dof_index: 3, index 1: 3, index 2: -1
+DEAL::  index 4 global_dof_index: 1, index 1: -1, index 2: 0
+DEAL::  index 5 global_dof_index: 4, index 1: -1, index 2: 1
+DEAL::  index 6 global_dof_index: 3, index 1: -1, index 2: 2
+DEAL::  index 7 global_dof_index: 5, index 1: -1, index 2: 3
+DEAL::n_quadrature_points: 16
+DEAL::  q index 0, index 1: 0, index 2: 0, q1: 0.211325 0.211325, q2: 1.21132 0.211325
+DEAL::  q index 1, index 1: 0, index 2: 1, q1: 0.211325 0.211325, q2: 1.78868 0.211325
+DEAL::  q index 2, index 1: 0, index 2: 2, q1: 0.211325 0.211325, q2: 1.21132 0.788675
+DEAL::  q index 3, index 1: 0, index 2: 3, q1: 0.211325 0.211325, q2: 1.78868 0.788675
+DEAL::  q index 4, index 1: 1, index 2: 0, q1: 0.788675 0.211325, q2: 1.21132 0.211325
+DEAL::  q index 5, index 1: 1, index 2: 1, q1: 0.788675 0.211325, q2: 1.78868 0.211325
+DEAL::  q index 6, index 1: 1, index 2: 2, q1: 0.788675 0.211325, q2: 1.21132 0.788675
+DEAL::  q index 7, index 1: 1, index 2: 3, q1: 0.788675 0.211325, q2: 1.78868 0.788675
+DEAL::  q index 8, index 1: 2, index 2: 0, q1: 0.211325 0.788675, q2: 1.21132 0.211325
+DEAL::  q index 9, index 1: 2, index 2: 1, q1: 0.211325 0.788675, q2: 1.78868 0.211325
+DEAL::  q index 10, index 1: 2, index 2: 2, q1: 0.211325 0.788675, q2: 1.21132 0.788675
+DEAL::  q index 11, index 1: 2, index 2: 3, q1: 0.211325 0.788675, q2: 1.78868 0.788675
+DEAL::  q index 12, index 1: 3, index 2: 0, q1: 0.788675 0.788675, q2: 1.21132 0.211325
+DEAL::  q index 13, index 1: 3, index 2: 1, q1: 0.788675 0.788675, q2: 1.78868 0.211325
+DEAL::  q index 14, index 1: 3, index 2: 2, q1: 0.788675 0.788675, q2: 1.21132 0.788675
+DEAL::  q index 15, index 1: 3, index 2: 3, q1: 0.788675 0.788675, q2: 1.78868 0.788675
+DEAL::
+DEAL::** coupling between cell1 and cell2 (unrolled) **
+DEAL::n_coupling_dofs: 8
+DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 
+DEAL::  index 0 global_dof_index: 0, index 1: 0, index 2: -1
+DEAL::  index 1 global_dof_index: 1, index 1: 1, index 2: -1
+DEAL::  index 2 global_dof_index: 2, index 1: 2, index 2: -1
+DEAL::  index 3 global_dof_index: 3, index 1: 3, index 2: -1
+DEAL::  index 4 global_dof_index: 1, index 1: -1, index 2: 0
+DEAL::  index 5 global_dof_index: 4, index 1: -1, index 2: 1
+DEAL::  index 6 global_dof_index: 3, index 1: -1, index 2: 2
+DEAL::  index 7 global_dof_index: 5, index 1: -1, index 2: 3
+DEAL::n_quadrature_points: 4
+DEAL::  q index 0, index 1: 0, index 2: 0, q1: 0.211325 0.211325, q2: 1.21132 0.211325
+DEAL::  q index 1, index 1: 1, index 2: 1, q1: 0.788675 0.211325, q2: 1.78868 0.211325
+DEAL::  q index 2, index 1: 2, index 2: 2, q1: 0.211325 0.788675, q2: 1.21132 0.788675
+DEAL::  q index 3, index 1: 3, index 2: 3, q1: 0.788675 0.788675, q2: 1.78868 0.788675
+DEAL::
+DEAL::** coupling between cell1 and cell2 Trapez (overlapping) **
+DEAL::n_coupling_dofs: 8
+DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 
+DEAL::  index 0 global_dof_index: 0, index 1: 0, index 2: -1
+DEAL::  index 1 global_dof_index: 1, index 1: 1, index 2: -1
+DEAL::  index 2 global_dof_index: 2, index 1: 2, index 2: -1
+DEAL::  index 3 global_dof_index: 3, index 1: 3, index 2: -1
+DEAL::  index 4 global_dof_index: 1, index 1: -1, index 2: 0
+DEAL::  index 5 global_dof_index: 4, index 1: -1, index 2: 1
+DEAL::  index 6 global_dof_index: 3, index 1: -1, index 2: 2
+DEAL::  index 7 global_dof_index: 5, index 1: -1, index 2: 3
+DEAL::n_quadrature_points: 2
+DEAL::  q index 0, index 1: 1, index 2: 0, q1: 1.00000 0.00000, q2: 1.00000 0.00000
+DEAL::  q index 1, index 1: 3, index 2: 2, q1: 1.00000 1.00000, q2: 1.00000 1.00000
+DEAL::
+DEAL::** coupling between cell1 and boundary face, trapez (overlapping) **
+DEAL::n_coupling_dofs: 8
+DEAL::coupling_dof_indices: 0 1 2 3 0 1 2 3 
+DEAL::  index 0 global_dof_index: 0, index 1: 0, index 2: -1
+DEAL::  index 1 global_dof_index: 1, index 1: 1, index 2: -1
+DEAL::  index 2 global_dof_index: 2, index 1: 2, index 2: -1
+DEAL::  index 3 global_dof_index: 3, index 1: 3, index 2: -1
+DEAL::  index 4 global_dof_index: 0, index 1: -1, index 2: 0
+DEAL::  index 5 global_dof_index: 1, index 1: -1, index 2: 1
+DEAL::  index 6 global_dof_index: 2, index 1: -1, index 2: 2
+DEAL::  index 7 global_dof_index: 3, index 1: -1, index 2: 3
+DEAL::n_quadrature_points: 2
+DEAL::  q index 0, index 1: 0, index 2: 0, q1: 0.00000 0.00000, q2: 0.00000 0.00000
+DEAL::  q index 1, index 1: 2, index 2: 1, q1: 0.00000 1.00000, q2: 0.00000 1.00000
+DEAL::
+DEAL::** coupling between cell1 and cell1->internal_face, trapez (overlapping) **
+DEAL::n_coupling_dofs: 8
+DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 
+DEAL::  index 0 global_dof_index: 0, index 1: 0, index 2: -1
+DEAL::  index 1 global_dof_index: 1, index 1: 1, index 2: -1
+DEAL::  index 2 global_dof_index: 2, index 1: 2, index 2: -1
+DEAL::  index 3 global_dof_index: 3, index 1: 3, index 2: -1
+DEAL::  index 4 global_dof_index: 1, index 1: -1, index 2: 0
+DEAL::  index 5 global_dof_index: 4, index 1: -1, index 2: 1
+DEAL::  index 6 global_dof_index: 3, index 1: -1, index 2: 2
+DEAL::  index 7 global_dof_index: 5, index 1: -1, index 2: 3
+DEAL::n_quadrature_points: 2
+DEAL::  q index 0, index 1: 1, index 2: 0, q1: 1.00000 0.00000, q2: 1.00000 0.00000
+DEAL::  q index 1, index 1: 3, index 2: 1, q1: 1.00000 1.00000, q2: 1.00000 1.00000
+DEAL::
+DEAL::** coupling between cell1 and cell2->internal_face, trapez (overlapping) **
+DEAL::n_coupling_dofs: 8
+DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 
+DEAL::  index 0 global_dof_index: 0, index 1: 0, index 2: -1
+DEAL::  index 1 global_dof_index: 1, index 1: 1, index 2: -1
+DEAL::  index 2 global_dof_index: 2, index 1: 2, index 2: -1
+DEAL::  index 3 global_dof_index: 3, index 1: 3, index 2: -1
+DEAL::  index 4 global_dof_index: 1, index 1: -1, index 2: 0
+DEAL::  index 5 global_dof_index: 4, index 1: -1, index 2: 1
+DEAL::  index 6 global_dof_index: 3, index 1: -1, index 2: 2
+DEAL::  index 7 global_dof_index: 5, index 1: -1, index 2: 3
+DEAL::n_quadrature_points: 2
+DEAL::  q index 0, index 1: 1, index 2: 0, q1: 1.00000 0.00000, q2: 1.00000 0.00000
+DEAL::  q index 1, index 1: 3, index 2: 1, q1: 1.00000 1.00000, q2: 1.00000 1.00000
+DEAL::
+DEAL::** coupling between cell1 and cell2->other_face, trapez (overlapping) **
+DEAL::n_coupling_dofs: 8
+DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 
+DEAL::  index 0 global_dof_index: 0, index 1: 0, index 2: -1
+DEAL::  index 1 global_dof_index: 1, index 1: 1, index 2: -1
+DEAL::  index 2 global_dof_index: 2, index 1: 2, index 2: -1
+DEAL::  index 3 global_dof_index: 3, index 1: 3, index 2: -1
+DEAL::  index 4 global_dof_index: 1, index 1: -1, index 2: 0
+DEAL::  index 5 global_dof_index: 4, index 1: -1, index 2: 1
+DEAL::  index 6 global_dof_index: 3, index 1: -1, index 2: 2
+DEAL::  index 7 global_dof_index: 5, index 1: -1, index 2: 3
+DEAL::n_quadrature_points: 0
+DEAL::
+DEAL::** coupling between cell1->face(0) and cell1->face(2), trapez (tensor) **
+DEAL::n_coupling_dofs: 8
+DEAL::coupling_dof_indices: 0 1 2 3 0 1 2 3 
+DEAL::  index 0 global_dof_index: 0, index 1: 0, index 2: -1
+DEAL::  index 1 global_dof_index: 1, index 1: 1, index 2: -1
+DEAL::  index 2 global_dof_index: 2, index 1: 2, index 2: -1
+DEAL::  index 3 global_dof_index: 3, index 1: 3, index 2: -1
+DEAL::  index 4 global_dof_index: 0, index 1: -1, index 2: 0
+DEAL::  index 5 global_dof_index: 1, index 1: -1, index 2: 1
+DEAL::  index 6 global_dof_index: 2, index 1: -1, index 2: 2
+DEAL::  index 7 global_dof_index: 3, index 1: -1, index 2: 3
+DEAL::n_quadrature_points: 4
+DEAL::  q index 0, index 1: 0, index 2: 0, q1: 0.00000 0.00000, q2: 0.00000 0.00000
+DEAL::  q index 1, index 1: 0, index 2: 1, q1: 0.00000 0.00000, q2: 1.00000 0.00000
+DEAL::  q index 2, index 1: 1, index 2: 0, q1: 0.00000 1.00000, q2: 0.00000 0.00000
+DEAL::  q index 3, index 1: 1, index 2: 1, q1: 0.00000 1.00000, q2: 1.00000 0.00000
+DEAL::
diff --git a/tests/fecoupling/fe_coupling_02.cc b/tests/fecoupling/fe_coupling_02.cc
new file mode 100644 (file)
index 0000000..3f589ac
--- /dev/null
@@ -0,0 +1,140 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test QuadratureCouplingType::tensor_product integration in FECouplingValues
+// with DoFCouplingType::independent. Typical for BEM coupling.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_coupling_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../test_grids.h"
+
+int
+main()
+{
+  initlog(0);
+  constexpr unsigned int dim = 1;
+
+  Triangulation<dim> tria;
+  TestGrids::hyper_line(tria, 2);
+
+  DoFHandler<dim> dofh(tria);
+  FE_Q<dim>       fe(1);
+  dofh.distribute_dofs(fe);
+
+  UpdateFlags update_flags =
+    update_quadrature_points | update_values | update_JxW_values;
+
+  FEValues<dim> fv1(fe, QGauss<dim>(fe.degree + 1), update_flags);
+  FEValues<dim> fv2(fe, QGauss<dim>(fe.degree + 1), update_flags);
+
+
+  std::vector<types::global_dof_index> v1(fe.dofs_per_cell);
+  std::vector<types::global_dof_index> v2(fe.dofs_per_cell);
+
+  const auto &cell1 = dofh.begin();
+  fv1.reinit(cell1);
+  cell1->get_dof_indices(v1);
+
+  unsigned int boundary_face = numbers::invalid_unsigned_int;
+  unsigned int internal_face = numbers::invalid_unsigned_int;
+
+  for (const unsigned int f : cell1->face_indices())
+    {
+      if (cell1->at_boundary(f))
+        {
+          boundary_face = f;
+        }
+      else if (!cell1->at_boundary(f))
+        {
+          internal_face = f;
+        }
+      if (boundary_face != numbers::invalid_unsigned_int &&
+          internal_face != numbers::invalid_unsigned_int)
+        break;
+    }
+
+  const auto &cell2 = cell1->neighbor(internal_face);
+  fv2.reinit(cell2);
+  cell2->get_dof_indices(v2);
+
+  // Now we have two cells, cell1 and cell2, and we want to couple them. We
+  // integrate over the tensor product of the two cells:
+  //
+  // \int_T1 \int_T2 |x1-x2| v_i(x1) v_j(x2) dx1 dx2
+  //
+
+  FECouplingValues<dim> fcv(fv1,
+                            fv2,
+                            DoFCouplingType::independent,
+                            QuadratureCouplingType::tensor_product);
+
+  deallog << "dofs_per_cell = " << fe.dofs_per_cell << std::endl
+          << "n_q_points_per_cell = " << fv2.n_quadrature_points << std::endl
+          << "n_q_points = " << fcv.n_quadrature_points() << std::endl;
+
+  FullMatrix<double> matrix(fcv.n_first_dofs(), fcv.n_second_dofs());
+
+  const auto first  = fcv.get_first_extractor(FEValuesExtractors::Scalar(0));
+  const auto second = fcv.get_second_extractor(FEValuesExtractors::Scalar(0));
+
+  // We need to loop over all the coupling quadrature points
+  for (const auto q : fcv.quadrature_point_indices())
+    {
+      const auto &[x, y] = fcv.quadrature_point(q);
+      for (const auto i : fcv.first_dof_indices())
+        {
+          const auto &vi = fcv[first].value(i, q);
+          for (const auto j : fcv.second_dof_indices())
+            {
+              const auto &vj = fcv[second].value(j, q);
+
+              deallog << std::left                                           //
+                      << "vi[" << i << "]: " << std::setw(10) << vi          //
+                      << "vj[" << j << "]: " << std::setw(10) << vj          //
+                      << "x[" << q << "]: " << std::setw(10) << x            //
+                      << "y[" << q << "]: " << std::setw(10) << y            //
+                      << "JxW[" << q << "]: " << std::setw(10) << fcv.JxW(q) //
+                      << "distance = " << x.distance(y) << std::endl;
+              matrix(i, j) += x.distance(y) * vi * vj * fcv.JxW(q);
+            }
+        }
+    }
+  // Now print the matrix:
+  deallog << "matrix = " << std::endl;
+  for (unsigned int i = 0; i < matrix.m(); ++i)
+    {
+      for (unsigned int j = 0; j < matrix.n(); ++j)
+        {
+          deallog << matrix(i, j) << " ";
+        }
+      deallog << std::endl;
+    }
+}
diff --git a/tests/fecoupling/fe_coupling_02.output b/tests/fecoupling/fe_coupling_02.output
new file mode 100644 (file)
index 0000000..bdddcd6
--- /dev/null
@@ -0,0 +1,23 @@
+
+DEAL::dofs_per_cell = 2
+DEAL::n_q_points_per_cell = 2
+DEAL::n_q_points = 4
+DEAL::vi[0]: 0.788675  vj[0]: 0.788675  x[0]: 0.211325  y[0]: 1.21132   JxW[0]: 0.250000  distance = 1.00000
+DEAL::vi[0]: 0.788675  vj[1]: 0.211325  x[0]: 0.211325  y[0]: 1.21132   JxW[0]: 0.250000  distance = 1.00000
+DEAL::vi[1]: 0.211325  vj[0]: 0.788675  x[0]: 0.211325  y[0]: 1.21132   JxW[0]: 0.250000  distance = 1.00000
+DEAL::vi[1]: 0.211325  vj[1]: 0.211325  x[0]: 0.211325  y[0]: 1.21132   JxW[0]: 0.250000  distance = 1.00000
+DEAL::vi[0]: 0.788675  vj[0]: 0.211325  x[1]: 0.211325  y[1]: 1.78868   JxW[1]: 0.250000  distance = 1.57735
+DEAL::vi[0]: 0.788675  vj[1]: 0.788675  x[1]: 0.211325  y[1]: 1.78868   JxW[1]: 0.250000  distance = 1.57735
+DEAL::vi[1]: 0.211325  vj[0]: 0.211325  x[1]: 0.211325  y[1]: 1.78868   JxW[1]: 0.250000  distance = 1.57735
+DEAL::vi[1]: 0.211325  vj[1]: 0.788675  x[1]: 0.211325  y[1]: 1.78868   JxW[1]: 0.250000  distance = 1.57735
+DEAL::vi[0]: 0.211325  vj[0]: 0.788675  x[2]: 0.788675  y[2]: 1.21132   JxW[2]: 0.250000  distance = 0.422650
+DEAL::vi[0]: 0.211325  vj[1]: 0.211325  x[2]: 0.788675  y[2]: 1.21132   JxW[2]: 0.250000  distance = 0.422650
+DEAL::vi[1]: 0.788675  vj[0]: 0.788675  x[2]: 0.788675  y[2]: 1.21132   JxW[2]: 0.250000  distance = 0.422650
+DEAL::vi[1]: 0.788675  vj[1]: 0.211325  x[2]: 0.788675  y[2]: 1.21132   JxW[2]: 0.250000  distance = 0.422650
+DEAL::vi[0]: 0.211325  vj[0]: 0.211325  x[3]: 0.788675  y[3]: 1.78868   JxW[3]: 0.250000  distance = 1.00000
+DEAL::vi[0]: 0.211325  vj[1]: 0.788675  x[3]: 0.788675  y[3]: 1.78868   JxW[3]: 0.250000  distance = 1.00000
+DEAL::vi[1]: 0.788675  vj[0]: 0.211325  x[3]: 0.788675  y[3]: 1.78868   JxW[3]: 0.250000  distance = 1.00000
+DEAL::vi[1]: 0.788675  vj[1]: 0.788675  x[3]: 0.788675  y[3]: 1.78868   JxW[3]: 0.250000  distance = 1.00000
+DEAL::matrix = 
+DEAL::0.250000 0.333333 
+DEAL::0.166667 0.250000 
diff --git a/tests/fecoupling/fe_coupling_03.cc b/tests/fecoupling/fe_coupling_03.cc
new file mode 100644 (file)
index 0000000..bdc4f2b
--- /dev/null
@@ -0,0 +1,169 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test QuadratureCouplingType::matching integration in FECouplingValues
+// with DoFCouplingType::contiguous. Typical for DG interfaces.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_coupling_values.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../test_grids.h"
+
+int
+main()
+{
+  initlog();
+  constexpr unsigned int dim = 2;
+
+  Triangulation<dim> tria;
+  TestGrids::hyper_line(tria, 2);
+
+  DoFHandler<dim> dofh(tria);
+  FE_DGQ<dim>     fe(1);
+  dofh.distribute_dofs(fe);
+
+  UpdateFlags update_flags =
+    update_quadrature_points | update_values | update_JxW_values;
+
+  QGauss<dim - 1> quad(fe.degree + 1);
+
+  FEFaceValues<dim>      fv1(fe, quad, update_flags);
+  FEFaceValues<dim>      fv2(fe, quad, update_flags);
+  FEInterfaceValues<dim> fiv(fe, quad, update_flags);
+
+  std::vector<types::global_dof_index> v1(fe.dofs_per_cell);
+  std::vector<types::global_dof_index> v2(fe.dofs_per_cell);
+
+  const auto &cell1 = dofh.begin();
+  cell1->get_dof_indices(v1);
+
+  unsigned int boundary_face = numbers::invalid_unsigned_int;
+  unsigned int internal_face = numbers::invalid_unsigned_int;
+
+  for (const unsigned int f : cell1->face_indices())
+    {
+      if (cell1->at_boundary(f))
+        {
+          boundary_face = f;
+        }
+      else if (!cell1->at_boundary(f))
+        {
+          internal_face = f;
+        }
+      if (boundary_face != numbers::invalid_unsigned_int &&
+          internal_face != numbers::invalid_unsigned_int)
+        break;
+    }
+  fv1.reinit(cell1, internal_face);
+
+  const auto &cell2 = cell1->neighbor(internal_face);
+  cell2->get_dof_indices(v2);
+  fv2.reinit(cell2, cell1->neighbor_of_neighbor(internal_face));
+
+  // Now reinit the interface values
+  fiv.reinit(cell1,
+             internal_face,
+             numbers::invalid_unsigned_int,
+             cell2,
+             cell1->neighbor_of_neighbor(internal_face),
+             numbers::invalid_unsigned_int);
+
+  // Now we have two cells, cell1 and cell2, and we want to integrate on the
+  // face.
+  //
+  // \int_f [v_i](x1) * [v_j](x2)) dx1 dx2
+  //
+
+  FECouplingValues<dim> fcv(fv1,
+                            fv2,
+                            DoFCouplingType::contiguous,
+                            QuadratureCouplingType::matching);
+
+  deallog << "dofs_per_cell = " << fe.dofs_per_cell << std::endl
+          << "n_q_points_per_cell = " << fv2.n_quadrature_points << std::endl
+          << "n_q_points = " << fcv.n_quadrature_points() << std::endl
+          << "n_interface_dofs = " << fiv.n_current_interface_dofs()
+          << std::endl
+          << "n_coupling_dofs = " << fcv.n_coupling_dofs() << std::endl;
+
+  FullMatrix<double> matrix(fcv.n_coupling_dofs(), fcv.n_coupling_dofs());
+  FEValuesExtractors::Scalar scalar(0);
+
+  const auto first  = fcv.get_first_extractor(scalar);
+  const auto second = fcv.get_second_extractor(scalar);
+
+  // We need to loop over all the coupling quadrature points
+  for (const auto q : fcv.quadrature_point_indices())
+    {
+      const auto &[x, y] = fcv.quadrature_point(q);
+      for (const auto i : fcv.coupling_dof_indices())
+        {
+          const auto &lvi       = fcv[first].value(i, q);
+          const auto &rvi       = fcv[second].value(i, q);
+          const auto  jump_i    = lvi - rvi;
+          const auto &jump_ifiv = fiv[scalar].jump_in_values(i, q);
+
+
+          for (const auto j : fcv.coupling_dof_indices())
+            {
+              const auto &lvj = fcv[first].value(j, q);
+              const auto &rvj = fcv[second].value(j, q);
+
+              const auto  jump_j    = lvj - rvj;
+              const auto &jump_jfiv = fiv[scalar].jump_in_values(j, q);
+
+              deallog << std::left << "fcv -- "                           //
+                      << "[vi][" << i << "]: " << std::setw(10) << jump_i //
+                      << "[vj][" << j << "]: " << std::setw(10) << jump_j //
+                      << "x[" << q << "]: " << std::setw(10) << x         //
+                      << "  y[" << q << "]: " << std::setw(10) << y       //
+                      << "JxW[" << q << "]: " << std::setw(10) << fcv.JxW(q)
+                      << std::endl;
+              deallog << std::left << "fiv -- "                              //
+                      << "[vi][" << i << "]: " << std::setw(10) << jump_ifiv //
+                      << "[vj][" << j << "]: " << std::setw(10) << jump_jfiv //
+                      << std::endl;
+              matrix(i, j) += jump_i * jump_j * fcv.JxW(q);
+            }
+        }
+    }
+  // Now print the matrix:
+  deallog << "matrix = " << std::endl;
+  for (unsigned int i = 0; i < matrix.m(); ++i)
+    {
+      for (unsigned int j = 0; j < matrix.n(); ++j)
+        {
+          deallog << matrix(i, j) << " ";
+        }
+      deallog << std::endl;
+    }
+}
diff --git a/tests/fecoupling/fe_coupling_03.output b/tests/fecoupling/fe_coupling_03.output
new file mode 100644 (file)
index 0000000..14ce313
--- /dev/null
@@ -0,0 +1,271 @@
+
+DEAL::dofs_per_cell = 4
+DEAL::n_q_points_per_cell = 2
+DEAL::n_q_points = 2
+DEAL::n_interface_dofs = 8
+DEAL::n_coupling_dofs = 8
+DEAL::fcv -- [vi][0]: -0.00000  [vj][0]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][0]: 0.00000   
+DEAL::fcv -- [vi][0]: -0.00000  [vj][1]: 0.788675  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][1]: 0.788675  
+DEAL::fcv -- [vi][0]: -0.00000  [vj][2]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][2]: 0.00000   
+DEAL::fcv -- [vi][0]: -0.00000  [vj][3]: 0.211325  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][3]: 0.211325  
+DEAL::fcv -- [vi][0]: -0.00000  [vj][4]: -0.788675 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][4]: -0.788675 
+DEAL::fcv -- [vi][0]: -0.00000  [vj][5]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][5]: 0.00000   
+DEAL::fcv -- [vi][0]: -0.00000  [vj][6]: -0.211325 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][6]: -0.211325 
+DEAL::fcv -- [vi][0]: -0.00000  [vj][7]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][7]: 0.00000   
+DEAL::fcv -- [vi][1]: 0.788675  [vj][0]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.788675  [vj][0]: 0.00000   
+DEAL::fcv -- [vi][1]: 0.788675  [vj][1]: 0.788675  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.788675  [vj][1]: 0.788675  
+DEAL::fcv -- [vi][1]: 0.788675  [vj][2]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.788675  [vj][2]: 0.00000   
+DEAL::fcv -- [vi][1]: 0.788675  [vj][3]: 0.211325  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.788675  [vj][3]: 0.211325  
+DEAL::fcv -- [vi][1]: 0.788675  [vj][4]: -0.788675 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.788675  [vj][4]: -0.788675 
+DEAL::fcv -- [vi][1]: 0.788675  [vj][5]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.788675  [vj][5]: 0.00000   
+DEAL::fcv -- [vi][1]: 0.788675  [vj][6]: -0.211325 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.788675  [vj][6]: -0.211325 
+DEAL::fcv -- [vi][1]: 0.788675  [vj][7]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.788675  [vj][7]: 0.00000   
+DEAL::fcv -- [vi][2]: -0.00000  [vj][0]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][0]: 0.00000   
+DEAL::fcv -- [vi][2]: -0.00000  [vj][1]: 0.788675  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][1]: 0.788675  
+DEAL::fcv -- [vi][2]: -0.00000  [vj][2]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][2]: 0.00000   
+DEAL::fcv -- [vi][2]: -0.00000  [vj][3]: 0.211325  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][3]: 0.211325  
+DEAL::fcv -- [vi][2]: -0.00000  [vj][4]: -0.788675 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][4]: -0.788675 
+DEAL::fcv -- [vi][2]: -0.00000  [vj][5]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][5]: 0.00000   
+DEAL::fcv -- [vi][2]: -0.00000  [vj][6]: -0.211325 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][6]: -0.211325 
+DEAL::fcv -- [vi][2]: -0.00000  [vj][7]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][7]: 0.00000   
+DEAL::fcv -- [vi][3]: 0.211325  [vj][0]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.211325  [vj][0]: 0.00000   
+DEAL::fcv -- [vi][3]: 0.211325  [vj][1]: 0.788675  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.211325  [vj][1]: 0.788675  
+DEAL::fcv -- [vi][3]: 0.211325  [vj][2]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.211325  [vj][2]: 0.00000   
+DEAL::fcv -- [vi][3]: 0.211325  [vj][3]: 0.211325  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.211325  [vj][3]: 0.211325  
+DEAL::fcv -- [vi][3]: 0.211325  [vj][4]: -0.788675 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.211325  [vj][4]: -0.788675 
+DEAL::fcv -- [vi][3]: 0.211325  [vj][5]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.211325  [vj][5]: 0.00000   
+DEAL::fcv -- [vi][3]: 0.211325  [vj][6]: -0.211325 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.211325  [vj][6]: -0.211325 
+DEAL::fcv -- [vi][3]: 0.211325  [vj][7]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.211325  [vj][7]: 0.00000   
+DEAL::fcv -- [vi][4]: -0.788675 [vj][0]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.788675 [vj][0]: 0.00000   
+DEAL::fcv -- [vi][4]: -0.788675 [vj][1]: 0.788675  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.788675 [vj][1]: 0.788675  
+DEAL::fcv -- [vi][4]: -0.788675 [vj][2]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.788675 [vj][2]: 0.00000   
+DEAL::fcv -- [vi][4]: -0.788675 [vj][3]: 0.211325  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.788675 [vj][3]: 0.211325  
+DEAL::fcv -- [vi][4]: -0.788675 [vj][4]: -0.788675 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.788675 [vj][4]: -0.788675 
+DEAL::fcv -- [vi][4]: -0.788675 [vj][5]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.788675 [vj][5]: 0.00000   
+DEAL::fcv -- [vi][4]: -0.788675 [vj][6]: -0.211325 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.788675 [vj][6]: -0.211325 
+DEAL::fcv -- [vi][4]: -0.788675 [vj][7]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.788675 [vj][7]: 0.00000   
+DEAL::fcv -- [vi][5]: 0.00000   [vj][0]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][0]: 0.00000   
+DEAL::fcv -- [vi][5]: 0.00000   [vj][1]: 0.788675  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][1]: 0.788675  
+DEAL::fcv -- [vi][5]: 0.00000   [vj][2]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][2]: 0.00000   
+DEAL::fcv -- [vi][5]: 0.00000   [vj][3]: 0.211325  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][3]: 0.211325  
+DEAL::fcv -- [vi][5]: 0.00000   [vj][4]: -0.788675 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][4]: -0.788675 
+DEAL::fcv -- [vi][5]: 0.00000   [vj][5]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][5]: 0.00000   
+DEAL::fcv -- [vi][5]: 0.00000   [vj][6]: -0.211325 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][6]: -0.211325 
+DEAL::fcv -- [vi][5]: 0.00000   [vj][7]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][7]: 0.00000   
+DEAL::fcv -- [vi][6]: -0.211325 [vj][0]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.211325 [vj][0]: 0.00000   
+DEAL::fcv -- [vi][6]: -0.211325 [vj][1]: 0.788675  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.211325 [vj][1]: 0.788675  
+DEAL::fcv -- [vi][6]: -0.211325 [vj][2]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.211325 [vj][2]: 0.00000   
+DEAL::fcv -- [vi][6]: -0.211325 [vj][3]: 0.211325  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.211325 [vj][3]: 0.211325  
+DEAL::fcv -- [vi][6]: -0.211325 [vj][4]: -0.788675 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.211325 [vj][4]: -0.788675 
+DEAL::fcv -- [vi][6]: -0.211325 [vj][5]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.211325 [vj][5]: 0.00000   
+DEAL::fcv -- [vi][6]: -0.211325 [vj][6]: -0.211325 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.211325 [vj][6]: -0.211325 
+DEAL::fcv -- [vi][6]: -0.211325 [vj][7]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.211325 [vj][7]: 0.00000   
+DEAL::fcv -- [vi][7]: 0.00000   [vj][0]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][0]: 0.00000   
+DEAL::fcv -- [vi][7]: 0.00000   [vj][1]: 0.788675  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][1]: 0.788675  
+DEAL::fcv -- [vi][7]: 0.00000   [vj][2]: -0.00000  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][2]: 0.00000   
+DEAL::fcv -- [vi][7]: 0.00000   [vj][3]: 0.211325  x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][3]: 0.211325  
+DEAL::fcv -- [vi][7]: 0.00000   [vj][4]: -0.788675 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][4]: -0.788675 
+DEAL::fcv -- [vi][7]: 0.00000   [vj][5]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][5]: 0.00000   
+DEAL::fcv -- [vi][7]: 0.00000   [vj][6]: -0.211325 x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][6]: -0.211325 
+DEAL::fcv -- [vi][7]: 0.00000   [vj][7]: 0.00000   x[0]: 1.00000    0.211325  y[0]: 1.00000    0.211325JxW[0]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][7]: 0.00000   
+DEAL::fcv -- [vi][0]: -0.00000  [vj][0]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][0]: 0.00000   
+DEAL::fcv -- [vi][0]: -0.00000  [vj][1]: 0.211325  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][1]: 0.211325  
+DEAL::fcv -- [vi][0]: -0.00000  [vj][2]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][2]: 0.00000   
+DEAL::fcv -- [vi][0]: -0.00000  [vj][3]: 0.788675  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][3]: 0.788675  
+DEAL::fcv -- [vi][0]: -0.00000  [vj][4]: -0.211325 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][4]: -0.211325 
+DEAL::fcv -- [vi][0]: -0.00000  [vj][5]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][5]: 0.00000   
+DEAL::fcv -- [vi][0]: -0.00000  [vj][6]: -0.788675 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][6]: -0.788675 
+DEAL::fcv -- [vi][0]: -0.00000  [vj][7]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][0]: 0.00000   [vj][7]: 0.00000   
+DEAL::fcv -- [vi][1]: 0.211325  [vj][0]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.211325  [vj][0]: 0.00000   
+DEAL::fcv -- [vi][1]: 0.211325  [vj][1]: 0.211325  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.211325  [vj][1]: 0.211325  
+DEAL::fcv -- [vi][1]: 0.211325  [vj][2]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.211325  [vj][2]: 0.00000   
+DEAL::fcv -- [vi][1]: 0.211325  [vj][3]: 0.788675  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.211325  [vj][3]: 0.788675  
+DEAL::fcv -- [vi][1]: 0.211325  [vj][4]: -0.211325 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.211325  [vj][4]: -0.211325 
+DEAL::fcv -- [vi][1]: 0.211325  [vj][5]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.211325  [vj][5]: 0.00000   
+DEAL::fcv -- [vi][1]: 0.211325  [vj][6]: -0.788675 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.211325  [vj][6]: -0.788675 
+DEAL::fcv -- [vi][1]: 0.211325  [vj][7]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][1]: 0.211325  [vj][7]: 0.00000   
+DEAL::fcv -- [vi][2]: -0.00000  [vj][0]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][0]: 0.00000   
+DEAL::fcv -- [vi][2]: -0.00000  [vj][1]: 0.211325  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][1]: 0.211325  
+DEAL::fcv -- [vi][2]: -0.00000  [vj][2]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][2]: 0.00000   
+DEAL::fcv -- [vi][2]: -0.00000  [vj][3]: 0.788675  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][3]: 0.788675  
+DEAL::fcv -- [vi][2]: -0.00000  [vj][4]: -0.211325 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][4]: -0.211325 
+DEAL::fcv -- [vi][2]: -0.00000  [vj][5]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][5]: 0.00000   
+DEAL::fcv -- [vi][2]: -0.00000  [vj][6]: -0.788675 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][6]: -0.788675 
+DEAL::fcv -- [vi][2]: -0.00000  [vj][7]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][2]: 0.00000   [vj][7]: 0.00000   
+DEAL::fcv -- [vi][3]: 0.788675  [vj][0]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.788675  [vj][0]: 0.00000   
+DEAL::fcv -- [vi][3]: 0.788675  [vj][1]: 0.211325  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.788675  [vj][1]: 0.211325  
+DEAL::fcv -- [vi][3]: 0.788675  [vj][2]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.788675  [vj][2]: 0.00000   
+DEAL::fcv -- [vi][3]: 0.788675  [vj][3]: 0.788675  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.788675  [vj][3]: 0.788675  
+DEAL::fcv -- [vi][3]: 0.788675  [vj][4]: -0.211325 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.788675  [vj][4]: -0.211325 
+DEAL::fcv -- [vi][3]: 0.788675  [vj][5]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.788675  [vj][5]: 0.00000   
+DEAL::fcv -- [vi][3]: 0.788675  [vj][6]: -0.788675 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.788675  [vj][6]: -0.788675 
+DEAL::fcv -- [vi][3]: 0.788675  [vj][7]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][3]: 0.788675  [vj][7]: 0.00000   
+DEAL::fcv -- [vi][4]: -0.211325 [vj][0]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.211325 [vj][0]: 0.00000   
+DEAL::fcv -- [vi][4]: -0.211325 [vj][1]: 0.211325  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.211325 [vj][1]: 0.211325  
+DEAL::fcv -- [vi][4]: -0.211325 [vj][2]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.211325 [vj][2]: 0.00000   
+DEAL::fcv -- [vi][4]: -0.211325 [vj][3]: 0.788675  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.211325 [vj][3]: 0.788675  
+DEAL::fcv -- [vi][4]: -0.211325 [vj][4]: -0.211325 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.211325 [vj][4]: -0.211325 
+DEAL::fcv -- [vi][4]: -0.211325 [vj][5]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.211325 [vj][5]: 0.00000   
+DEAL::fcv -- [vi][4]: -0.211325 [vj][6]: -0.788675 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.211325 [vj][6]: -0.788675 
+DEAL::fcv -- [vi][4]: -0.211325 [vj][7]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][4]: -0.211325 [vj][7]: 0.00000   
+DEAL::fcv -- [vi][5]: 0.00000   [vj][0]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][0]: 0.00000   
+DEAL::fcv -- [vi][5]: 0.00000   [vj][1]: 0.211325  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][1]: 0.211325  
+DEAL::fcv -- [vi][5]: 0.00000   [vj][2]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][2]: 0.00000   
+DEAL::fcv -- [vi][5]: 0.00000   [vj][3]: 0.788675  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][3]: 0.788675  
+DEAL::fcv -- [vi][5]: 0.00000   [vj][4]: -0.211325 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][4]: -0.211325 
+DEAL::fcv -- [vi][5]: 0.00000   [vj][5]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][5]: 0.00000   
+DEAL::fcv -- [vi][5]: 0.00000   [vj][6]: -0.788675 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][6]: -0.788675 
+DEAL::fcv -- [vi][5]: 0.00000   [vj][7]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][5]: 0.00000   [vj][7]: 0.00000   
+DEAL::fcv -- [vi][6]: -0.788675 [vj][0]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.788675 [vj][0]: 0.00000   
+DEAL::fcv -- [vi][6]: -0.788675 [vj][1]: 0.211325  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.788675 [vj][1]: 0.211325  
+DEAL::fcv -- [vi][6]: -0.788675 [vj][2]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.788675 [vj][2]: 0.00000   
+DEAL::fcv -- [vi][6]: -0.788675 [vj][3]: 0.788675  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.788675 [vj][3]: 0.788675  
+DEAL::fcv -- [vi][6]: -0.788675 [vj][4]: -0.211325 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.788675 [vj][4]: -0.211325 
+DEAL::fcv -- [vi][6]: -0.788675 [vj][5]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.788675 [vj][5]: 0.00000   
+DEAL::fcv -- [vi][6]: -0.788675 [vj][6]: -0.788675 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.788675 [vj][6]: -0.788675 
+DEAL::fcv -- [vi][6]: -0.788675 [vj][7]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][6]: -0.788675 [vj][7]: 0.00000   
+DEAL::fcv -- [vi][7]: 0.00000   [vj][0]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][0]: 0.00000   
+DEAL::fcv -- [vi][7]: 0.00000   [vj][1]: 0.211325  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][1]: 0.211325  
+DEAL::fcv -- [vi][7]: 0.00000   [vj][2]: -0.00000  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][2]: 0.00000   
+DEAL::fcv -- [vi][7]: 0.00000   [vj][3]: 0.788675  x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][3]: 0.788675  
+DEAL::fcv -- [vi][7]: 0.00000   [vj][4]: -0.211325 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][4]: -0.211325 
+DEAL::fcv -- [vi][7]: 0.00000   [vj][5]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][5]: 0.00000   
+DEAL::fcv -- [vi][7]: 0.00000   [vj][6]: -0.788675 x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][6]: -0.788675 
+DEAL::fcv -- [vi][7]: 0.00000   [vj][7]: 0.00000   x[1]: 1.00000    0.788675  y[1]: 1.00000    0.788675JxW[1]: 0.500000  
+DEAL::fiv -- [vi][7]: 0.00000   [vj][7]: 0.00000   
+DEAL::matrix = 
+DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::0.00000 0.333333 0.00000 0.166667 -0.333333 0.00000 -0.166667 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::0.00000 0.166667 0.00000 0.333333 -0.166667 0.00000 -0.333333 0.00000 
+DEAL::0.00000 -0.333333 0.00000 -0.166667 0.333333 0.00000 0.166667 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
+DEAL::0.00000 -0.166667 0.00000 -0.333333 0.166667 0.00000 0.333333 0.00000 
+DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 
diff --git a/tests/fecoupling/fe_coupling_04.cc b/tests/fecoupling/fe_coupling_04.cc
new file mode 100644 (file)
index 0000000..2342960
--- /dev/null
@@ -0,0 +1,154 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test QuadratureCouplingType::reorder integration in FECouplingValues
+// with DoFCouplingType::contiguous. Typical for bulk-surface coupling.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_coupling_values.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../test_grids.h"
+
+int
+main()
+{
+  initlog();
+  constexpr unsigned int dim   = 2;
+  constexpr unsigned int codim = dim - 1;
+
+  Triangulation<dim> tria;
+  TestGrids::hyper_line(tria, 2);
+
+  Triangulation<codim, dim> codim_tria;
+  TestGrids::hyper_line(codim_tria, 2);
+
+  DoFHandler<dim> dofh(tria);
+  FE_Q<dim>       fe(1);
+  dofh.distribute_dofs(fe);
+
+  DoFHandler<codim, dim> codim_dofh(codim_tria);
+  FE_Q<codim, dim>       codim_fe(1);
+  codim_dofh.distribute_dofs(codim_fe);
+
+  UpdateFlags update_flags =
+    update_quadrature_points | update_values | update_JxW_values;
+
+  QGauss<codim> quad(fe.degree + 1);
+
+  FEFaceValues<dim>    fv1(fe, quad, update_flags);
+  FEValues<codim, dim> fv2(codim_fe, quad, update_flags);
+
+  std::vector<types::global_dof_index> v1(fe.dofs_per_cell);
+  std::vector<types::global_dof_index> v2(codim_fe.dofs_per_cell);
+
+  deallog << "Bulk dofs: " << fe.dofs_per_cell
+          << ", Surface dofs: " << codim_fe.dofs_per_cell << std::endl;
+
+  const auto &cell1 = dofh.begin();
+  cell1->get_dof_indices(v1);
+  // Coincides with the first cell of codim tria.
+  const unsigned int face = 2;
+
+  const auto &cell2 = codim_dofh.begin();
+  cell2->get_dof_indices(v2);
+
+  fv1.reinit(cell1, face);
+  fv2.reinit(cell2);
+
+  // Now we have two cells, cell1 and cell2, and we want to integrate on the
+  // face.
+  //
+  // \int_f v_i(x) * v_j(x) dx
+  //
+
+  FECouplingValues<dim, codim, dim> fcv(fv1,
+                                        fv2,
+                                        DoFCouplingType::contiguous,
+                                        QuadratureCouplingType::reorder);
+
+  FullMatrix<double> matrix(fcv.n_coupling_dofs(), fcv.n_coupling_dofs());
+  FEValuesExtractors::Scalar scalar(0);
+
+  const auto bulk    = fcv.get_first_extractor(scalar);
+  const auto surface = fcv.get_second_extractor(scalar);
+
+  deallog << "Renumbering vectors: " << std::endl;
+
+  for (const auto i : fcv.coupling_dof_indices())
+    {
+      auto id_first  = fcv.coupling_dof_to_dof_indices(i).first;
+      auto id_second = fcv.coupling_dof_to_dof_indices(i).second;
+      deallog << "i: " << i << " bulk: " << (int)id_first
+              << " surface: " << (int)id_second << std::endl;
+    }
+
+  // We need to loop over all the coupling quadrature points
+  for (const auto q : fcv.quadrature_point_indices())
+    {
+      const auto &[x, y] = fcv.quadrature_point(q);
+      for (const auto i : fcv.coupling_dof_indices())
+        {
+          const auto &bulk_vi    = fcv[bulk].value(i, q);
+          const auto &surface_vi = fcv[surface].value(i, q);
+
+          for (const auto j : fcv.coupling_dof_indices())
+            {
+              const auto &bulk_vj    = fcv[bulk].value(j, q);
+              const auto &surface_vj = fcv[surface].value(j, q);
+
+              deallog << std::left                                           //
+                      << "bvi[" << i << "]: " << std::setw(10) << bulk_vi    //
+                      << "bvj[" << j << "]: " << std::setw(10) << bulk_vj    //
+                      << "svi[" << i << "]: " << std::setw(10) << surface_vi //
+                      << "svj[" << j << "]: " << std::setw(10) << surface_vj //
+                      << "x[" << q << "]: " << std::setw(10) << x            //
+                      << "  y[" << q << "]: " << std::setw(10) << y          //
+                      << std::endl;
+              matrix(i, j) += (bulk_vi * bulk_vj +             //
+                               10 * bulk_vi * surface_vj +     //
+                               100 * surface_vi * bulk_vj +    //
+                               1000 * surface_vi * surface_vj) //
+                              * fcv.JxW(q);
+            }
+        }
+    }
+  // Now print the matrix:
+  deallog << "matrix = " << std::endl << std::left;
+  for (unsigned int i = 0; i < matrix.m(); ++i)
+    {
+      for (unsigned int j = 0; j < matrix.n(); ++j)
+        {
+          deallog << std::setw(8) << matrix(i, j) << " ";
+        }
+      deallog << std::endl;
+    }
+}
diff --git a/tests/fecoupling/fe_coupling_04.output b/tests/fecoupling/fe_coupling_04.output
new file mode 100644 (file)
index 0000000..b182b41
--- /dev/null
@@ -0,0 +1,88 @@
+
+DEAL::Bulk dofs: 4, Surface dofs: 2
+DEAL::Renumbering vectors: 
+DEAL::i: 0 bulk: 0 surface: -1
+DEAL::i: 1 bulk: 1 surface: -1
+DEAL::i: 2 bulk: 2 surface: -1
+DEAL::i: 3 bulk: 3 surface: -1
+DEAL::i: 4 bulk: -1 surface: 0
+DEAL::i: 5 bulk: -1 surface: 1
+DEAL::bvi[0]: 0.788675  bvj[0]: 0.788675  svi[0]: 0.00000   svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.788675  bvj[1]: 0.211325  svi[0]: 0.00000   svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.788675  bvj[2]: 0.00000   svi[0]: 0.00000   svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.788675  bvj[3]: 0.00000   svi[0]: 0.00000   svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.788675  bvj[4]: 0.00000   svi[0]: 0.00000   svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.788675  bvj[5]: 0.00000   svi[0]: 0.00000   svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[0]: 0.788675  svi[1]: 0.00000   svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[1]: 0.211325  svi[1]: 0.00000   svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[2]: 0.00000   svi[1]: 0.00000   svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[3]: 0.00000   svi[1]: 0.00000   svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[4]: 0.00000   svi[1]: 0.00000   svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[5]: 0.00000   svi[1]: 0.00000   svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[0]: 0.788675  svi[2]: 0.00000   svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[1]: 0.211325  svi[2]: 0.00000   svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[2]: 0.00000   svi[2]: 0.00000   svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[3]: 0.00000   svi[2]: 0.00000   svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[4]: 0.00000   svi[2]: 0.00000   svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[5]: 0.00000   svi[2]: 0.00000   svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[0]: 0.788675  svi[3]: 0.00000   svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[1]: 0.211325  svi[3]: 0.00000   svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[2]: 0.00000   svi[3]: 0.00000   svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[3]: 0.00000   svi[3]: 0.00000   svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[4]: 0.00000   svi[3]: 0.00000   svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[5]: 0.00000   svi[3]: 0.00000   svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[0]: 0.788675  svi[4]: 0.788675  svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[1]: 0.211325  svi[4]: 0.788675  svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[2]: 0.00000   svi[4]: 0.788675  svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[3]: 0.00000   svi[4]: 0.788675  svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[4]: 0.00000   svi[4]: 0.788675  svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[5]: 0.00000   svi[4]: 0.788675  svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[0]: 0.788675  svi[5]: 0.211325  svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[1]: 0.211325  svi[5]: 0.211325  svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[2]: 0.00000   svi[5]: 0.211325  svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[3]: 0.00000   svi[5]: 0.211325  svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[4]: 0.00000   svi[5]: 0.211325  svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[5]: 0.00000   svi[5]: 0.211325  svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.211325  bvj[0]: 0.211325  svi[0]: 0.00000   svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[0]: 0.211325  bvj[1]: 0.788675  svi[0]: 0.00000   svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[0]: 0.211325  bvj[2]: 0.00000   svi[0]: 0.00000   svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[0]: 0.211325  bvj[3]: 0.00000   svi[0]: 0.00000   svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[0]: 0.211325  bvj[4]: 0.00000   svi[0]: 0.00000   svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[0]: 0.211325  bvj[5]: 0.00000   svi[0]: 0.00000   svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[0]: 0.211325  svi[1]: 0.00000   svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[1]: 0.788675  svi[1]: 0.00000   svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[2]: 0.00000   svi[1]: 0.00000   svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[3]: 0.00000   svi[1]: 0.00000   svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[4]: 0.00000   svi[1]: 0.00000   svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[5]: 0.00000   svi[1]: 0.00000   svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[0]: 0.211325  svi[2]: 0.00000   svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[1]: 0.788675  svi[2]: 0.00000   svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[2]: 0.00000   svi[2]: 0.00000   svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[3]: 0.00000   svi[2]: 0.00000   svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[4]: 0.00000   svi[2]: 0.00000   svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[5]: 0.00000   svi[2]: 0.00000   svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[0]: 0.211325  svi[3]: 0.00000   svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[1]: 0.788675  svi[3]: 0.00000   svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[2]: 0.00000   svi[3]: 0.00000   svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[3]: 0.00000   svi[3]: 0.00000   svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[4]: 0.00000   svi[3]: 0.00000   svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[5]: 0.00000   svi[3]: 0.00000   svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[0]: 0.211325  svi[4]: 0.211325  svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[1]: 0.788675  svi[4]: 0.211325  svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[2]: 0.00000   svi[4]: 0.211325  svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[3]: 0.00000   svi[4]: 0.211325  svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[4]: 0.00000   svi[4]: 0.211325  svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[5]: 0.00000   svi[4]: 0.211325  svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[0]: 0.211325  svi[5]: 0.788675  svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[1]: 0.788675  svi[5]: 0.788675  svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[2]: 0.00000   svi[5]: 0.788675  svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[3]: 0.00000   svi[5]: 0.788675  svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[4]: 0.00000   svi[5]: 0.788675  svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[5]: 0.00000   svi[5]: 0.788675  svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::matrix = 
+DEAL::0.333333 0.166667 0.00000  0.00000  3.33333  1.66667  
+DEAL::0.166667 0.333333 0.00000  0.00000  1.66667  3.33333  
+DEAL::0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  
+DEAL::0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  
+DEAL::33.3333  16.6667  0.00000  0.00000  333.333  166.667  
+DEAL::16.6667  33.3333  0.00000  0.00000  166.667  333.333  
diff --git a/tests/fecoupling/fe_coupling_05.cc b/tests/fecoupling/fe_coupling_05.cc
new file mode 100644 (file)
index 0000000..b83e428
--- /dev/null
@@ -0,0 +1,154 @@
+// ---------------------------------------------------------------------
+//
+// Copyright (C) 2023 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.
+//
+// ---------------------------------------------------------------------
+
+
+// Test QuadratureCouplingType::matching integration in FECouplingValues
+// with DoFCouplingType::contiguous. Typical for bulk-surface coupling.
+
+#include <deal.II/base/quadrature_lib.h>
+
+#include <deal.II/fe/fe_coupling_values.h>
+#include <deal.II/fe/fe_dgq.h>
+#include <deal.II/fe/fe_interface_values.h>
+#include <deal.II/fe/fe_q.h>
+#include <deal.II/fe/fe_values.h>
+#include <deal.II/fe/mapping_q.h>
+
+#include <deal.II/grid/grid_generator.h>
+#include <deal.II/grid/grid_refinement.h>
+#include <deal.II/grid/tria.h>
+
+#include <fstream>
+#include <iomanip>
+#include <iostream>
+
+#include "../tests.h"
+
+#include "../test_grids.h"
+
+int
+main()
+{
+  initlog();
+  constexpr unsigned int dim   = 2;
+  constexpr unsigned int codim = dim - 1;
+
+  Triangulation<dim> tria;
+  TestGrids::hyper_line(tria, 2);
+
+  Triangulation<codim, dim> codim_tria;
+  TestGrids::hyper_line(codim_tria, 2);
+
+  DoFHandler<dim> dofh(tria);
+  FE_Q<dim>       fe(1);
+  dofh.distribute_dofs(fe);
+
+  DoFHandler<codim, dim> codim_dofh(codim_tria);
+  FE_Q<codim, dim>       codim_fe(1);
+  codim_dofh.distribute_dofs(codim_fe);
+
+  UpdateFlags update_flags =
+    update_quadrature_points | update_values | update_JxW_values;
+
+  QGauss<codim> quad(fe.degree + 1);
+
+  FEFaceValues<dim>    fv1(fe, quad, update_flags);
+  FEValues<codim, dim> fv2(codim_fe, quad, update_flags);
+
+  std::vector<types::global_dof_index> v1(fe.dofs_per_cell);
+  std::vector<types::global_dof_index> v2(codim_fe.dofs_per_cell);
+
+  deallog << "Bulk dofs: " << fe.dofs_per_cell
+          << ", Surface dofs: " << codim_fe.dofs_per_cell << std::endl;
+
+  const auto &cell1 = dofh.begin();
+  cell1->get_dof_indices(v1);
+  // Coincides with the first cell of codim tria.
+  const unsigned int face = 2;
+
+  const auto &cell2 = codim_dofh.begin();
+  cell2->get_dof_indices(v2);
+
+  fv1.reinit(cell1, face);
+  fv2.reinit(cell2);
+
+  // Now we have two cells, cell1 and cell2, and we want to integrate on the
+  // face.
+  //
+  // \int_f v_i(x) * v_j(x) dx
+  //
+
+  FECouplingValues<dim, codim, dim> fcv(fv1,
+                                        fv2,
+                                        DoFCouplingType::contiguous,
+                                        QuadratureCouplingType::matching);
+
+  FullMatrix<double> matrix(fcv.n_coupling_dofs(), fcv.n_coupling_dofs());
+  FEValuesExtractors::Scalar scalar(0);
+
+  const auto bulk    = fcv.get_first_extractor(scalar);
+  const auto surface = fcv.get_second_extractor(scalar);
+
+  deallog << "Renumbering vectors: " << std::endl;
+
+  for (const auto i : fcv.coupling_dof_indices())
+    {
+      auto id_first  = fcv.coupling_dof_to_dof_indices(i).first;
+      auto id_second = fcv.coupling_dof_to_dof_indices(i).second;
+      deallog << "i: " << i << " bulk: " << (int)id_first
+              << " surface: " << (int)id_second << std::endl;
+    }
+
+  // We need to loop over all the coupling quadrature points
+  for (const auto q : fcv.quadrature_point_indices())
+    {
+      const auto &[x, y] = fcv.quadrature_point(q);
+      for (const auto i : fcv.coupling_dof_indices())
+        {
+          const auto &bulk_vi    = fcv[bulk].value(i, q);
+          const auto &surface_vi = fcv[surface].value(i, q);
+
+          for (const auto j : fcv.coupling_dof_indices())
+            {
+              const auto &bulk_vj    = fcv[bulk].value(j, q);
+              const auto &surface_vj = fcv[surface].value(j, q);
+
+              deallog << std::left                                           //
+                      << "bvi[" << i << "]: " << std::setw(10) << bulk_vi    //
+                      << "bvj[" << j << "]: " << std::setw(10) << bulk_vj    //
+                      << "svi[" << i << "]: " << std::setw(10) << surface_vi //
+                      << "svj[" << j << "]: " << std::setw(10) << surface_vj //
+                      << "x[" << q << "]: " << std::setw(10) << x            //
+                      << "  y[" << q << "]: " << std::setw(10) << y          //
+                      << std::endl;
+              matrix(i, j) += (bulk_vi * bulk_vj +             //
+                               10 * bulk_vi * surface_vj +     //
+                               100 * surface_vi * bulk_vj +    //
+                               1000 * surface_vi * surface_vj) //
+                              * fcv.JxW(q);
+            }
+        }
+    }
+  // Now print the matrix:
+  deallog << "matrix = " << std::endl << std::left;
+  for (unsigned int i = 0; i < matrix.m(); ++i)
+    {
+      for (unsigned int j = 0; j < matrix.n(); ++j)
+        {
+          deallog << std::setw(8) << matrix(i, j) << " ";
+        }
+      deallog << std::endl;
+    }
+}
diff --git a/tests/fecoupling/fe_coupling_05.output b/tests/fecoupling/fe_coupling_05.output
new file mode 100644 (file)
index 0000000..b182b41
--- /dev/null
@@ -0,0 +1,88 @@
+
+DEAL::Bulk dofs: 4, Surface dofs: 2
+DEAL::Renumbering vectors: 
+DEAL::i: 0 bulk: 0 surface: -1
+DEAL::i: 1 bulk: 1 surface: -1
+DEAL::i: 2 bulk: 2 surface: -1
+DEAL::i: 3 bulk: 3 surface: -1
+DEAL::i: 4 bulk: -1 surface: 0
+DEAL::i: 5 bulk: -1 surface: 1
+DEAL::bvi[0]: 0.788675  bvj[0]: 0.788675  svi[0]: 0.00000   svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.788675  bvj[1]: 0.211325  svi[0]: 0.00000   svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.788675  bvj[2]: 0.00000   svi[0]: 0.00000   svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.788675  bvj[3]: 0.00000   svi[0]: 0.00000   svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.788675  bvj[4]: 0.00000   svi[0]: 0.00000   svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.788675  bvj[5]: 0.00000   svi[0]: 0.00000   svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[0]: 0.788675  svi[1]: 0.00000   svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[1]: 0.211325  svi[1]: 0.00000   svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[2]: 0.00000   svi[1]: 0.00000   svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[3]: 0.00000   svi[1]: 0.00000   svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[4]: 0.00000   svi[1]: 0.00000   svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[1]: 0.211325  bvj[5]: 0.00000   svi[1]: 0.00000   svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[0]: 0.788675  svi[2]: 0.00000   svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[1]: 0.211325  svi[2]: 0.00000   svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[2]: 0.00000   svi[2]: 0.00000   svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[3]: 0.00000   svi[2]: 0.00000   svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[4]: 0.00000   svi[2]: 0.00000   svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[2]: 0.00000   bvj[5]: 0.00000   svi[2]: 0.00000   svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[0]: 0.788675  svi[3]: 0.00000   svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[1]: 0.211325  svi[3]: 0.00000   svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[2]: 0.00000   svi[3]: 0.00000   svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[3]: 0.00000   svi[3]: 0.00000   svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[4]: 0.00000   svi[3]: 0.00000   svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[3]: 0.00000   bvj[5]: 0.00000   svi[3]: 0.00000   svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[0]: 0.788675  svi[4]: 0.788675  svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[1]: 0.211325  svi[4]: 0.788675  svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[2]: 0.00000   svi[4]: 0.788675  svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[3]: 0.00000   svi[4]: 0.788675  svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[4]: 0.00000   svi[4]: 0.788675  svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[4]: 0.00000   bvj[5]: 0.00000   svi[4]: 0.788675  svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[0]: 0.788675  svi[5]: 0.211325  svj[0]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[1]: 0.211325  svi[5]: 0.211325  svj[1]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[2]: 0.00000   svi[5]: 0.211325  svj[2]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[3]: 0.00000   svi[5]: 0.211325  svj[3]: 0.00000   x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[4]: 0.00000   svi[5]: 0.211325  svj[4]: 0.788675  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[5]: 0.00000   bvj[5]: 0.00000   svi[5]: 0.211325  svj[5]: 0.211325  x[0]: 0.211325   0.00000  y[0]: 0.211325   0.00000
+DEAL::bvi[0]: 0.211325  bvj[0]: 0.211325  svi[0]: 0.00000   svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[0]: 0.211325  bvj[1]: 0.788675  svi[0]: 0.00000   svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[0]: 0.211325  bvj[2]: 0.00000   svi[0]: 0.00000   svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[0]: 0.211325  bvj[3]: 0.00000   svi[0]: 0.00000   svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[0]: 0.211325  bvj[4]: 0.00000   svi[0]: 0.00000   svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[0]: 0.211325  bvj[5]: 0.00000   svi[0]: 0.00000   svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[0]: 0.211325  svi[1]: 0.00000   svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[1]: 0.788675  svi[1]: 0.00000   svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[2]: 0.00000   svi[1]: 0.00000   svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[3]: 0.00000   svi[1]: 0.00000   svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[4]: 0.00000   svi[1]: 0.00000   svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[1]: 0.788675  bvj[5]: 0.00000   svi[1]: 0.00000   svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[0]: 0.211325  svi[2]: 0.00000   svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[1]: 0.788675  svi[2]: 0.00000   svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[2]: 0.00000   svi[2]: 0.00000   svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[3]: 0.00000   svi[2]: 0.00000   svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[4]: 0.00000   svi[2]: 0.00000   svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[2]: 0.00000   bvj[5]: 0.00000   svi[2]: 0.00000   svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[0]: 0.211325  svi[3]: 0.00000   svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[1]: 0.788675  svi[3]: 0.00000   svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[2]: 0.00000   svi[3]: 0.00000   svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[3]: 0.00000   svi[3]: 0.00000   svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[4]: 0.00000   svi[3]: 0.00000   svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[3]: 0.00000   bvj[5]: 0.00000   svi[3]: 0.00000   svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[0]: 0.211325  svi[4]: 0.211325  svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[1]: 0.788675  svi[4]: 0.211325  svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[2]: 0.00000   svi[4]: 0.211325  svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[3]: 0.00000   svi[4]: 0.211325  svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[4]: 0.00000   svi[4]: 0.211325  svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[4]: 0.00000   bvj[5]: 0.00000   svi[4]: 0.211325  svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[0]: 0.211325  svi[5]: 0.788675  svj[0]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[1]: 0.788675  svi[5]: 0.788675  svj[1]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[2]: 0.00000   svi[5]: 0.788675  svj[2]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[3]: 0.00000   svi[5]: 0.788675  svj[3]: 0.00000   x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[4]: 0.00000   svi[5]: 0.788675  svj[4]: 0.211325  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::bvi[5]: 0.00000   bvj[5]: 0.00000   svi[5]: 0.788675  svj[5]: 0.788675  x[1]: 0.788675   0.00000  y[1]: 0.788675   0.00000
+DEAL::matrix = 
+DEAL::0.333333 0.166667 0.00000  0.00000  3.33333  1.66667  
+DEAL::0.166667 0.333333 0.00000  0.00000  1.66667  3.33333  
+DEAL::0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  
+DEAL::0.00000  0.00000  0.00000  0.00000  0.00000  0.00000  
+DEAL::33.3333  16.6667  0.00000  0.00000  333.333  166.667  
+DEAL::16.6667  33.3333  0.00000  0.00000  166.667  333.333  

In the beginning the Universe was created. This has made a lot of people very angry and has been widely regarded as a bad move.

Douglas Adams


Typeset in Trocchi and Trocchi Bold Sans Serif.