--- /dev/null
+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)
#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>
} // 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
}
} // 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
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
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
--- /dev/null
+cmake_minimum_required(VERSION 3.13.4)
+include(../scripts/setup_testsubproject.cmake)
+project(testsuite CXX)
+deal_ii_pickup_tests()
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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>();
+}
--- /dev/null
+
+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::
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ }
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ }
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ }
+}
--- /dev/null
+
+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
--- /dev/null
+// ---------------------------------------------------------------------
+//
+// 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;
+ }
+}
--- /dev/null
+
+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
projects / dealii.git / commitdiff