From 588ac3dbcec3f8b19f84f82d856d172a39af87e7 Mon Sep 17 00:00:00 2001 From: Luca Heltai Date: Fri, 21 Jul 2023 00:54:48 +0200 Subject: [PATCH] FECouplingValues. --- doc/news/changes/minor/20230721LucaHeltai | 4 + include/deal.II/fe/fe_coupling_values.h | 1158 +++++++++++++++++++++ include/deal.II/fe/fe_values_extractors.h | 49 +- tests/fecoupling/CMakeLists.txt | 4 + tests/fecoupling/fe_coupling_01.cc | 224 ++++ tests/fecoupling/fe_coupling_01.output | 137 +++ tests/fecoupling/fe_coupling_02.cc | 140 +++ tests/fecoupling/fe_coupling_02.output | 23 + tests/fecoupling/fe_coupling_03.cc | 169 +++ tests/fecoupling/fe_coupling_03.output | 271 +++++ tests/fecoupling/fe_coupling_04.cc | 154 +++ tests/fecoupling/fe_coupling_04.output | 88 ++ tests/fecoupling/fe_coupling_05.cc | 154 +++ tests/fecoupling/fe_coupling_05.output | 88 ++ 14 files changed, 2662 insertions(+), 1 deletion(-) create mode 100644 doc/news/changes/minor/20230721LucaHeltai create mode 100644 tests/fecoupling/CMakeLists.txt create mode 100644 tests/fecoupling/fe_coupling_01.cc create mode 100644 tests/fecoupling/fe_coupling_01.output create mode 100644 tests/fecoupling/fe_coupling_02.cc create mode 100644 tests/fecoupling/fe_coupling_02.output create mode 100644 tests/fecoupling/fe_coupling_03.cc create mode 100644 tests/fecoupling/fe_coupling_03.output create mode 100644 tests/fecoupling/fe_coupling_04.cc create mode 100644 tests/fecoupling/fe_coupling_04.output create mode 100644 tests/fecoupling/fe_coupling_05.cc create mode 100644 tests/fecoupling/fe_coupling_05.output diff --git a/doc/news/changes/minor/20230721LucaHeltai b/doc/news/changes/minor/20230721LucaHeltai new file mode 100644 index 0000000000..9caf84be9c --- /dev/null +++ b/doc/news/changes/minor/20230721LucaHeltai @@ -0,0 +1,4 @@ +New: Created a new class FECouplingValues, that helps with assembling of coupled +finite element methods across different dimensions or different grids. +
+(Luca Heltai, 2023/07/21) diff --git a/include/deal.II/fe/fe_coupling_values.h b/include/deal.II/fe/fe_coupling_values.h index c6ab804fec..2ed108beec 100644 --- a/include/deal.II/fe/fe_coupling_values.h +++ b/include/deal.II/fe/fe_coupling_values.h @@ -20,6 +20,7 @@ #include +#include #include #include #include @@ -406,6 +407,664 @@ namespace FEValuesViews } // namespace FEValuesViews +/** + * Quadrature coupling options when assembling quadrature formulas for double + * integrals. + * + * When computing the approximation of double integrals of the form + * + * \f[ + * \int_{T_1} \int{T_2} K(x_1, x_2) f(x_1) g(x_2) dT_1 dT_2, + * \f] + * + * where $T_1$ and $T_2$ are two arbitrary sets (cells, faces, edges, or any + * combination thereof), and $K$ is a (possibly singular) coupling kernel, one + * needs to combine quadrature formulas from two different FEValuesBase objects. + * + * This enum class provides a way to specify how the quadrature points and + * weights should be combined. In general, the two FEValuesBase objects provide + * different quadrature rules, and these can be interpreted in different ways, + * depending on the kernel function that is being integrated, and on how the two + * quadrature rules were constructed. + * + * This enum is used in the constructor of FECouplingValues to specify how to + * interpret and manipulate the quadrature points and weights of the two + * FEValuesBase objects. + */ +enum class QuadratureCouplingType +{ + /** + * The FEValuesBase objects provide different quadrature rules, and the + * resulting quadrature points and weights should be constructed as the tensor + * product of the two quadrature rules. This is generally used for non-local + * and non-singular kernels, where the quadrature points of the two + * FEValuesBase objects are independent of each other. + */ + tensor_product, + + /** + * Both FEValuesBase objects provide the same number of (generally different) + * quadrature points, that should be used as is to implement the double + * integration. This is equivalent to rewriting double integrals as a single + * sum over unrolled indices, and it is useful, for example, when performing + * integration of singular kernels, which require special quadrature rules + * both on the inside and on the outside integral. An exception is thrown if + * the two FEValuesBase objects do not provide the same number of quadrature + * points, but otherwise the two sets of points and weights can be arbitrary. + */ + unrolled, + + /** + * The FEValuesBase objects provide the same quadrature rule over the same + * set, with the same orientation. This is similar to the unrolled case, in + * the sense that the quadrature formulas are used as is, but both objects are + * assumed to return the same quadrature points and weights. This is useful + * when integrating over faces of neighboring cells, and when the reordering + * of the quadrature points is known to match a priori. An exception is thrown + * if the quadrature points and weights of the two FEValuesBase objects do not + * match one-to-one. + */ + matching, + + /** + * The FEValuesBase objects provide the same quadrature rule over the same + * set, but with possibly different orientations. The quadrature points are + * reordered internally, so that the resulting quadrature points and weights + * are the same as in the matching case. An exception is thrown if the + * quadrature points and weights of the two FEValuesBase cannot be reordered + * to match one-to-one. + */ + reorder, + + /** + * The FEValuesBase objects provide partially overlapping quadrature rules + * over two intersecting sets, with possibly different orientations. The + * quadrature points are reordered and matched internally, so that the + * resulting quadrature points and weights are only the ones that match + * between the two objects. If no overlapping points and weights can be found, + * an empty set is used to integrate, and the resulting integral is zero. No + * exception is thrown in this case. + */ + overlapping, +}; + +/** + * DoF coupling options when assembling double integrals. + * + * When computing the approximation of double integrals of the form + * + * \f[ + * \int_{T_1} \int{T_2} K(x_1, x_2) v_i(x_1) w_j(x_2) dT_1 dT_2, + * \f] + * + * where $T_1$ and $T_2$ are two arbitrary sets (cells, faces, edges, or any + * combination thereof), and $K$ is a (possibly singular) coupling kernel, one + * may want to combine degrees from two different FEValuesBase objects (i.e., + * basis functions $v_i$ and $w_j$ in the examples above). + * + * This enum class provides a way to specify how the degrees of freedom should + * be combined. There are two cases of interest: + * + * 1. the two FEValuesBase objects refer to different DoFHandlers + * 2. the two FEValuesBase objects refer to the same DoFHandler + * + * In the first case, one usually treats the two sets of degrees of freedom as + * independent of each other, and the resulting matrix is generally rectangular. + * + * In the second case, one may choose to treat the two sets of degrees of + * freedom either as independent or to group them together. A similar approach + * is used in the FEInterfaceValues class, where the degrees of freedom of the + * two FEValuesBase objects are grouped together, in a contiguous way, so that + * the resulting basis functions are interpreted in the following way: + * + * \f[ + * \phi_{1,i}(x) = \begin{cases} v_i(x) & \text{ if } i \in [0,n_l) \\ + * 0 & \text{ if ) i \in [n_1, n_1+n_2] \end{cases},\quad \phi_{1,i}(x) = + * \begin{cases} 0(x) & \text{ if } i \in [0,n_1) \\ + * w_{i-n_1}(x) & \text{ if ) i \in [n_1, n_1+n_2] \end{cases}, + * \f] + * + * where $phi_{1,i}$ is the first basis function with index $i$ and $n_{1,2}$ + * are the number of local dofs on the first and second FEValuesBase objects. + * + * This enum is used in the constructor of FECouplingValues to specify how to + * interpret and manipulate the local dof indices of the two FEValuesBase + * objects. + */ +enum class DoFCouplingType +{ + /** + * The FEValuesBase objects may have different dof indices, possibly indexing + * different DoFHandler objects, and we are interested in assembling a + * generally rectangular matrix, where there is no relationship between the + * two index spaces. + */ + independent, + + /** + * The FEValuesBase objects may have different dof indices, but they both + * index the same DoFHandler objects, and we are interested in assembling them + * all together in a single matrix. In this coupling type, the DoF indices are + * grouped together, one after the other, first the first dof indices, and + * then the second dof indices. This is useful when one wants to assemble four + * matrices at the same time, corresponding to the four possible combinations + * of the two FEValuesBase objects, (i.e., first coupled with first, first + * coupled with second, second coupled with first, and second coupled with + * second) and then sum them together to obtain the final system. This is + * similar to what is done in the FEInterfaceValues class. + */ + contiguous, +}; + +/** + * FECouplingValues is a class that facilitates the integration of finite + * element data between two different finite element objects, possibly living on + * different grids, and with possibly different topological dimensions (i.e., + * cells, faces, edges, and any combination thereof). + * + * This class provides a way to simplify the implementation of the following + * abstract operation: + * + * \f[ + * \int_{T_1} \int{T_2} K(x_1, x_2) \phi^1_i(x_1) \phi^2_j(x_2) dT_1 dT_2 + * \f] + * + * for three different types of Kernels $K$: + * - $K(x_1, x_2)$ is a non-singular Kernel function, for example, it is a + * function of positive powers $\alpha$ of the distance between the quadrature + * points $|x_1-x_2|^\alpha$; + * - $K(x_1, x_2)$ is a singular Kernel function, for example, it is a function + * of negative powers $\alpha$ of the distance between the quadrature points + * $|x_1-x_2|^\alpha$; + * - $K(x_1, x_2)$ is a Dirac delta distribution $\delta(x_1-x_2)$, such that + * the integral above is actually a single integral over the intersection of + * the two sets $T_1$ and $T_2$. + * + * For the first case, one may think that the only natural way to proceed is to + * compute the double integral by simply nesting two loops: + * \f[ + * \int_{T_1} \int{T_2} K(x_1, x_2) \phi^1_i(x_1) \phi^2_j(x_2) dT_1 dT_2 + * \approx \sum_{q_1} \sum_{q_2} K(x_1^{q_1}, x_2^{q_2}) \phi^1_i(x_1^{q_1}) + * \phi^2_j(x_2^{q_2}) w_1^{q_1} w_2^{q_2}, + * \f] + * + * where $x_1^{q_1}$ and $x_2^{q_2}$ are the quadrature points in $T_1$ and + * $T_2$ respectively, and $w_1^{q_1}$ and $w_2^{q_2}$ are the corresponding + * quadrature weights. + * + * This, however is not the only way to proceed. In fact, such an integral can + * be rewritten as a single loop over corresponding elements of two arrays of + * points with the same length that can be thought of as a single quadrature + * rule on the set $T_1\times T_2$. For singular kernels, for example, this is + * often the only way to proceed, since the quadrature formula on $T_1\times + * T_2$ is usually not written as a tensor product quadrature formula, and one + * needs to build a custom quadrature formula for this purpose. + * + * This class allows one to treat the three cases above in the same way, and to + * approximate the integral as follows: + * + * \f[ + * \int_{T_1} \int{T_2} K(x_1, x_2) \phi^1_i(x_1) \phi^2_j(x_2) dT_1 dT_2 + * \approx \sum_{i=1}^{N_q} K(x_1^{i}, x_2^{i}) \phi^1_i(x_1^{i}) + * \phi^2_j(x_2^{i}) w_1^{i} w_2^i, + * \f] + * + * Since the triple of objects $(\{q\}, \{w\}, \{\phi\})$ is usually provided by + * a class derived from the FEValuesBase class, this is the type that the class + * needs at construction time. $T_1$ and $T_2$ can be two arbitrary cells, + * faces, or edges belonging to possibly different meshes (or to meshes with + * different topological dimensions), $\phi^1_i$ and $\phi^2_j$ are basis + * functions defined on $T_1$ and $T_2$, respectively. + * + * The case of the Dirac distribution is when $T_1$ and $T_2$ + * correspond to the common face of two neighboring cells. In this case, this + * class provides a functionality which is similar to the FEInterfaceValues + * class, and gives you a way to access values of basis functions on the + * neighboring cells, as well as their gradients and Hessians, in a unified + * fashion, on the face. + * + * Similarly, this class can be used to couple bulk and surface meshes across + * the faces of the bulk mesh. In this case, the two FEValuesBase objects will + * have different topological dimension (i.e., one will be a cell in a + * co-dimension one triangulation, and the other a face of a bulk grid with + * co-dimension zero), and the QuadratureCouplingType argument is usually chosen + * to be QuadratureCouplingType::reorder, since the quadrature points of the two + * different FEValuesBase objects are not necessarily generated with the same + * ordering. + * + * The type of integral to compute is controlled by the QuadratureCouplingType + * argument (see the documentation of that enum class for more details), while + * the type degrees of freedom coupling is controlled by the DoFCouplingType + * argument (see the documentation of that enum class for more details). + * + * As an example usage of this class, consider the a bilinear form of the form: + * \f[ + * \int_{T_1} \int{T_2} K_1(x_1, x_2) v_i(x_1) u_j(x_2) dT_1 dT_2 + + * \int_{T_1} \int{T_2} K_2(x_1, x_2) p_i(x_1) q_j(x_2) dT_1 dT_2 + * \f] + * where the finite dimensional space has two scalar components. We indicate + * with $v_i$ and $p_i$ the trial functions, and with $u_j$ and + * $q_j$ the corresponding test functions. $K_1$ and $K_2$ are coupling kernels: + * such a formulation is used, for example, to write the bilinear forms of + * Galerkin boundary element methods. + * + * The corresponding implementation would look like the following: + * + * @code + * ... // double loop over cells that yields cell_1 and cell_2 + * + * fe_values_1.reinit(cell_1); + * fe_values_2.reinit(cell_2); + * + * CouplingFEValues cfv(fe_values1, fe_values2, + * DoFCouplingType::independent, + * QuadratureCouplingType::tensor_product); + * + * FullMatrix 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 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 +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 &fe_values_1, + const FEValuesBase &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 &fe_values_1, + const FEValuesBase &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 + FEValuesExtractors::FirstCoupling + get_first_extractor(const Extractor &extractor) const; + + /** + * Return a SecondCoupling object that can be used to extract values from the + * second FEValuesBase object. + */ + template + FEValuesExtractors::SecondCoupling + 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> + 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 + 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 + 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 + 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 + 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 + get_coupling_dof_indices( + const std::vector &dof_indices_1, + const std::vector &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 + 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 + 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 + const FEValuesViews::RenumberedView< + FEValuesViews::View> + operator[]( + const FEValuesExtractors::FirstCoupling &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 + const FEValuesViews::RenumberedView< + FEValuesViews::View> + operator[]( + const FEValuesExtractors::SecondCoupling &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> first_fe_values; + + /** + * Pointer to second FEValuesBase object. + */ + SmartPointer> second_fe_values; + + /** + * Renumbering data for the first FEValuesBase object. + */ + std::unique_ptr first_renumbering_data; + + /** + * Renumbering data for the second FEValuesBase object. + */ + std::unique_ptr second_renumbering_data; + + /** + * Number of quadrature points. + */ + unsigned int n_quadrature_points_; + + /** + * Number of coupling DoF indices. If DoFCouplingType::independent is used, + * this is numbers::invalid_unsigned_int, while it is n_first_dofs() + + * n_second_dofs() in the the case of DoFCouplingType::contiguous. + */ + unsigned int n_coupling_dofs_; +}; + + #ifndef DOXYGEN @@ -660,6 +1319,505 @@ namespace FEValuesViews } } // namespace FEValuesViews +/*-------------- Inline functions FECouplingValues ---------------------*/ + +template +std::vector +FECouplingValues::get_coupling_dof_indices( + const std::vector &dof_indices_1, + const std::vector &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 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 +FECouplingValues::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 +FECouplingValues::FECouplingValues( + const FEValuesBase &fe_values_1, + const FEValuesBase &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 +void +FECouplingValues::reinit( + const FEValuesBase &fe_values_1, + const FEValuesBase &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 first_dofs_map; + std::vector second_dofs_map; + + std::vector first_quad_map; + std::vector 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( + first_n_inner_dofs, + first_n_inner_quadrature_points, + first_dofs_map, + first_quad_map); + + second_renumbering_data = std::make_unique( + second_n_inner_dofs, + second_n_inner_quadrature_points, + second_dofs_map, + second_quad_map); +} + + +template +inline double +FECouplingValues::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 +inline std_cxx20::ranges::iota_view +FECouplingValues::quadrature_point_indices() const +{ + return {0U, n_quadrature_points_}; +} + + + +template +inline std_cxx20::ranges::iota_view +FECouplingValues::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 +inline std_cxx20::ranges::iota_view +FECouplingValues::first_dof_indices() const +{ + return {0U, n_first_dofs()}; +} + + + +template +inline std_cxx20::ranges::iota_view +FECouplingValues::second_dof_indices() const +{ + return {0U, n_second_dofs()}; +} + + + +template +inline unsigned int +FECouplingValues::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 +inline unsigned int +FECouplingValues::n_first_dofs() const +{ + return first_renumbering_data->n_dofs; +} + + + +template +inline unsigned int +FECouplingValues::n_second_dofs() const +{ + return second_renumbering_data->n_dofs; +} + + + +template +inline unsigned int +FECouplingValues::n_quadrature_points() const +{ + return n_quadrature_points_; +} + + + +template +std::pair, Point> +FECouplingValues::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 +std::pair +FECouplingValues:: + 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 +std::pair +FECouplingValues::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 +template +inline FEValuesExtractors::FirstCoupling +FECouplingValues::get_first_extractor( + const Extractor &extractor) const +{ + return FEValuesExtractors::FirstCoupling(extractor); +} + + + +template +template +inline FEValuesExtractors::SecondCoupling +FECouplingValues::get_second_extractor( + const Extractor &extractor) const +{ + return FEValuesExtractors::SecondCoupling(extractor); +} + + + +template +template +inline const FEValuesViews::RenumberedView< + FEValuesViews::View> +FECouplingValues::operator[]( + const FEValuesExtractors::FirstCoupling &extractor) const +{ + return FEValuesViews::RenumberedView< + typename FEValuesViews::View>( + (*first_fe_values)[extractor.extractor], *first_renumbering_data); +} + + + +template +template +inline const FEValuesViews::RenumberedView< + typename FEValuesViews::View> +FECouplingValues::operator[]( + const FEValuesExtractors::SecondCoupling &extractor) const +{ + return FEValuesViews::RenumberedView< + typename FEValuesViews::View>( + (*second_fe_values)[extractor.extractor], *second_renumbering_data); +} + #endif // DOXYGEN DEAL_II_NAMESPACE_CLOSE diff --git a/include/deal.II/fe/fe_values_extractors.h b/include/deal.II/fe/fe_values_extractors.h index 2f73c16eb5..ef109bca01 100644 --- a/include/deal.II/fe/fe_values_extractors.h +++ b/include/deal.II/fe/fe_values_extractors.h @@ -268,6 +268,42 @@ namespace FEValuesExtractors std::string get_name() const; }; + + /** + * Helper struct to associate an extractor to the first FEValuesBase + * sub-object of a FECouplingValues object. + */ + template + 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 + struct SecondCoupling + { + /** + * Construct a new SecondCoupling object with the given extractor. + */ + SecondCoupling(const Extractor &extractor); + + /** + * The actual extractor object. + */ + const Extractor extractor; + }; } // namespace FEValuesExtractors @@ -320,9 +356,20 @@ namespace FEValuesExtractors inline Tensor::Tensor(const unsigned int first_tensor_component) : first_tensor_component(first_tensor_component) {} -} // namespace FEValuesExtractors + template + inline FirstCoupling::FirstCoupling(const Extractor &extractor) + : extractor(extractor) + {} + + + template + inline SecondCoupling::SecondCoupling(const Extractor &extractor) + : extractor(extractor) + {} +} // namespace FEValuesExtractors + DEAL_II_NAMESPACE_CLOSE diff --git a/tests/fecoupling/CMakeLists.txt b/tests/fecoupling/CMakeLists.txt new file mode 100644 index 0000000000..db9c1adcf2 --- /dev/null +++ b/tests/fecoupling/CMakeLists.txt @@ -0,0 +1,4 @@ +cmake_minimum_required(VERSION 3.13.4) +include(../scripts/setup_testsubproject.cmake) +project(testsuite CXX) +deal_ii_pickup_tests() diff --git a/tests/fecoupling/fe_coupling_01.cc b/tests/fecoupling/fe_coupling_01.cc new file mode 100644 index 0000000000..8ba61b76ce --- /dev/null +++ b/tests/fecoupling/fe_coupling_01.cc @@ -0,0 +1,224 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test basic properties of FECouplingValues + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include + +#include "../tests.h" + +#include "../test_grids.h" + + +template +void +inspect_fcv(const FCV &fcv, + const std::vector &v1, + const std::vector &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(pair.first) + << ", index 2: " << static_cast(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 +void +test() +{ + Triangulation tria; + TestGrids::hyper_line(tria, 2); + + DoFHandler dofh(tria); + FE_Q fe(1); + dofh.distribute_dofs(fe); + + MappingQ mapping(1); + UpdateFlags update_flags = update_quadrature_points | update_JxW_values; + + FEValues fv1(mapping, fe, QGauss(fe.degree + 1), update_flags); + FEValues fv2(mapping, fe, QGauss(fe.degree + 1), update_flags); + + FEValues fv3(mapping, fe, QTrapezoid(), update_flags); + FEValues fv4(mapping, fe, QTrapezoid(), update_flags); + + FEFaceValues fv5(mapping, fe, QTrapezoid(), update_flags); + FEFaceValues fv6(mapping, fe, QTrapezoid(), update_flags); + + std::vector v1(fe.dofs_per_cell); + std::vector 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 fcv(fv1, fv2, DoFCouplingType::contiguous); + inspect_fcv(fcv, v1, v2); + } + { + deallog << "** coupling between cell1 and cell2 (unrolled) **" << std::endl; + FECouplingValues 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 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 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 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 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 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 fcv(fv5, + fv6, + DoFCouplingType::contiguous, + QuadratureCouplingType::tensor_product); + inspect_fcv(fcv, v1, v1); + } +} + +int +main() +{ + initlog(); + test<2>(); +} diff --git a/tests/fecoupling/fe_coupling_01.output b/tests/fecoupling/fe_coupling_01.output new file mode 100644 index 0000000000..1e45d5224b --- /dev/null +++ b/tests/fecoupling/fe_coupling_01.output @@ -0,0 +1,137 @@ + +DEAL::** coupling between cell1 and cell2 (tensor) ** +DEAL::n_coupling_dofs: 8 +DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 +DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1 +DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1 +DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1 +DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1 +DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0 +DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1 +DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2 +DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3 +DEAL::n_quadrature_points: 16 +DEAL:: q index 0, index 1: 0, index 2: 0, q1: 0.211325 0.211325, q2: 1.21132 0.211325 +DEAL:: q index 1, index 1: 0, index 2: 1, q1: 0.211325 0.211325, q2: 1.78868 0.211325 +DEAL:: q index 2, index 1: 0, index 2: 2, q1: 0.211325 0.211325, q2: 1.21132 0.788675 +DEAL:: q index 3, index 1: 0, index 2: 3, q1: 0.211325 0.211325, q2: 1.78868 0.788675 +DEAL:: q index 4, index 1: 1, index 2: 0, q1: 0.788675 0.211325, q2: 1.21132 0.211325 +DEAL:: q index 5, index 1: 1, index 2: 1, q1: 0.788675 0.211325, q2: 1.78868 0.211325 +DEAL:: q index 6, index 1: 1, index 2: 2, q1: 0.788675 0.211325, q2: 1.21132 0.788675 +DEAL:: q index 7, index 1: 1, index 2: 3, q1: 0.788675 0.211325, q2: 1.78868 0.788675 +DEAL:: q index 8, index 1: 2, index 2: 0, q1: 0.211325 0.788675, q2: 1.21132 0.211325 +DEAL:: q index 9, index 1: 2, index 2: 1, q1: 0.211325 0.788675, q2: 1.78868 0.211325 +DEAL:: q index 10, index 1: 2, index 2: 2, q1: 0.211325 0.788675, q2: 1.21132 0.788675 +DEAL:: q index 11, index 1: 2, index 2: 3, q1: 0.211325 0.788675, q2: 1.78868 0.788675 +DEAL:: q index 12, index 1: 3, index 2: 0, q1: 0.788675 0.788675, q2: 1.21132 0.211325 +DEAL:: q index 13, index 1: 3, index 2: 1, q1: 0.788675 0.788675, q2: 1.78868 0.211325 +DEAL:: q index 14, index 1: 3, index 2: 2, q1: 0.788675 0.788675, q2: 1.21132 0.788675 +DEAL:: q index 15, index 1: 3, index 2: 3, q1: 0.788675 0.788675, q2: 1.78868 0.788675 +DEAL:: +DEAL::** coupling between cell1 and cell2 (unrolled) ** +DEAL::n_coupling_dofs: 8 +DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 +DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1 +DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1 +DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1 +DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1 +DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0 +DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1 +DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2 +DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3 +DEAL::n_quadrature_points: 4 +DEAL:: q index 0, index 1: 0, index 2: 0, q1: 0.211325 0.211325, q2: 1.21132 0.211325 +DEAL:: q index 1, index 1: 1, index 2: 1, q1: 0.788675 0.211325, q2: 1.78868 0.211325 +DEAL:: q index 2, index 1: 2, index 2: 2, q1: 0.211325 0.788675, q2: 1.21132 0.788675 +DEAL:: q index 3, index 1: 3, index 2: 3, q1: 0.788675 0.788675, q2: 1.78868 0.788675 +DEAL:: +DEAL::** coupling between cell1 and cell2 Trapez (overlapping) ** +DEAL::n_coupling_dofs: 8 +DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 +DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1 +DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1 +DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1 +DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1 +DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0 +DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1 +DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2 +DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3 +DEAL::n_quadrature_points: 2 +DEAL:: q index 0, index 1: 1, index 2: 0, q1: 1.00000 0.00000, q2: 1.00000 0.00000 +DEAL:: q index 1, index 1: 3, index 2: 2, q1: 1.00000 1.00000, q2: 1.00000 1.00000 +DEAL:: +DEAL::** coupling between cell1 and boundary face, trapez (overlapping) ** +DEAL::n_coupling_dofs: 8 +DEAL::coupling_dof_indices: 0 1 2 3 0 1 2 3 +DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1 +DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1 +DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1 +DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1 +DEAL:: index 4 global_dof_index: 0, index 1: -1, index 2: 0 +DEAL:: index 5 global_dof_index: 1, index 1: -1, index 2: 1 +DEAL:: index 6 global_dof_index: 2, index 1: -1, index 2: 2 +DEAL:: index 7 global_dof_index: 3, index 1: -1, index 2: 3 +DEAL::n_quadrature_points: 2 +DEAL:: q index 0, index 1: 0, index 2: 0, q1: 0.00000 0.00000, q2: 0.00000 0.00000 +DEAL:: q index 1, index 1: 2, index 2: 1, q1: 0.00000 1.00000, q2: 0.00000 1.00000 +DEAL:: +DEAL::** coupling between cell1 and cell1->internal_face, trapez (overlapping) ** +DEAL::n_coupling_dofs: 8 +DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 +DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1 +DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1 +DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1 +DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1 +DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0 +DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1 +DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2 +DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3 +DEAL::n_quadrature_points: 2 +DEAL:: q index 0, index 1: 1, index 2: 0, q1: 1.00000 0.00000, q2: 1.00000 0.00000 +DEAL:: q index 1, index 1: 3, index 2: 1, q1: 1.00000 1.00000, q2: 1.00000 1.00000 +DEAL:: +DEAL::** coupling between cell1 and cell2->internal_face, trapez (overlapping) ** +DEAL::n_coupling_dofs: 8 +DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 +DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1 +DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1 +DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1 +DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1 +DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0 +DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1 +DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2 +DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3 +DEAL::n_quadrature_points: 2 +DEAL:: q index 0, index 1: 1, index 2: 0, q1: 1.00000 0.00000, q2: 1.00000 0.00000 +DEAL:: q index 1, index 1: 3, index 2: 1, q1: 1.00000 1.00000, q2: 1.00000 1.00000 +DEAL:: +DEAL::** coupling between cell1 and cell2->other_face, trapez (overlapping) ** +DEAL::n_coupling_dofs: 8 +DEAL::coupling_dof_indices: 0 1 2 3 1 4 3 5 +DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1 +DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1 +DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1 +DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1 +DEAL:: index 4 global_dof_index: 1, index 1: -1, index 2: 0 +DEAL:: index 5 global_dof_index: 4, index 1: -1, index 2: 1 +DEAL:: index 6 global_dof_index: 3, index 1: -1, index 2: 2 +DEAL:: index 7 global_dof_index: 5, index 1: -1, index 2: 3 +DEAL::n_quadrature_points: 0 +DEAL:: +DEAL::** coupling between cell1->face(0) and cell1->face(2), trapez (tensor) ** +DEAL::n_coupling_dofs: 8 +DEAL::coupling_dof_indices: 0 1 2 3 0 1 2 3 +DEAL:: index 0 global_dof_index: 0, index 1: 0, index 2: -1 +DEAL:: index 1 global_dof_index: 1, index 1: 1, index 2: -1 +DEAL:: index 2 global_dof_index: 2, index 1: 2, index 2: -1 +DEAL:: index 3 global_dof_index: 3, index 1: 3, index 2: -1 +DEAL:: index 4 global_dof_index: 0, index 1: -1, index 2: 0 +DEAL:: index 5 global_dof_index: 1, index 1: -1, index 2: 1 +DEAL:: index 6 global_dof_index: 2, index 1: -1, index 2: 2 +DEAL:: index 7 global_dof_index: 3, index 1: -1, index 2: 3 +DEAL::n_quadrature_points: 4 +DEAL:: q index 0, index 1: 0, index 2: 0, q1: 0.00000 0.00000, q2: 0.00000 0.00000 +DEAL:: q index 1, index 1: 0, index 2: 1, q1: 0.00000 0.00000, q2: 1.00000 0.00000 +DEAL:: q index 2, index 1: 1, index 2: 0, q1: 0.00000 1.00000, q2: 0.00000 0.00000 +DEAL:: q index 3, index 1: 1, index 2: 1, q1: 0.00000 1.00000, q2: 1.00000 0.00000 +DEAL:: diff --git a/tests/fecoupling/fe_coupling_02.cc b/tests/fecoupling/fe_coupling_02.cc new file mode 100644 index 0000000000..3f589ac736 --- /dev/null +++ b/tests/fecoupling/fe_coupling_02.cc @@ -0,0 +1,140 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test QuadratureCouplingType::tensor_product integration in FECouplingValues +// with DoFCouplingType::independent. Typical for BEM coupling. + +#include + +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + +#include "../test_grids.h" + +int +main() +{ + initlog(0); + constexpr unsigned int dim = 1; + + Triangulation tria; + TestGrids::hyper_line(tria, 2); + + DoFHandler dofh(tria); + FE_Q fe(1); + dofh.distribute_dofs(fe); + + UpdateFlags update_flags = + update_quadrature_points | update_values | update_JxW_values; + + FEValues fv1(fe, QGauss(fe.degree + 1), update_flags); + FEValues fv2(fe, QGauss(fe.degree + 1), update_flags); + + + std::vector v1(fe.dofs_per_cell); + std::vector 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 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 matrix(fcv.n_first_dofs(), fcv.n_second_dofs()); + + const auto first = fcv.get_first_extractor(FEValuesExtractors::Scalar(0)); + const auto second = fcv.get_second_extractor(FEValuesExtractors::Scalar(0)); + + // We need to loop over all the coupling quadrature points + for (const auto q : fcv.quadrature_point_indices()) + { + const auto &[x, y] = fcv.quadrature_point(q); + for (const auto i : fcv.first_dof_indices()) + { + const auto &vi = fcv[first].value(i, q); + for (const auto j : fcv.second_dof_indices()) + { + const auto &vj = fcv[second].value(j, q); + + deallog << std::left // + << "vi[" << i << "]: " << std::setw(10) << vi // + << "vj[" << j << "]: " << std::setw(10) << vj // + << "x[" << q << "]: " << std::setw(10) << x // + << "y[" << q << "]: " << std::setw(10) << y // + << "JxW[" << q << "]: " << std::setw(10) << fcv.JxW(q) // + << "distance = " << x.distance(y) << std::endl; + matrix(i, j) += x.distance(y) * vi * vj * fcv.JxW(q); + } + } + } + // Now print the matrix: + deallog << "matrix = " << std::endl; + for (unsigned int i = 0; i < matrix.m(); ++i) + { + for (unsigned int j = 0; j < matrix.n(); ++j) + { + deallog << matrix(i, j) << " "; + } + deallog << std::endl; + } +} diff --git a/tests/fecoupling/fe_coupling_02.output b/tests/fecoupling/fe_coupling_02.output new file mode 100644 index 0000000000..bdddcd62a1 --- /dev/null +++ b/tests/fecoupling/fe_coupling_02.output @@ -0,0 +1,23 @@ + +DEAL::dofs_per_cell = 2 +DEAL::n_q_points_per_cell = 2 +DEAL::n_q_points = 4 +DEAL::vi[0]: 0.788675 vj[0]: 0.788675 x[0]: 0.211325 y[0]: 1.21132 JxW[0]: 0.250000 distance = 1.00000 +DEAL::vi[0]: 0.788675 vj[1]: 0.211325 x[0]: 0.211325 y[0]: 1.21132 JxW[0]: 0.250000 distance = 1.00000 +DEAL::vi[1]: 0.211325 vj[0]: 0.788675 x[0]: 0.211325 y[0]: 1.21132 JxW[0]: 0.250000 distance = 1.00000 +DEAL::vi[1]: 0.211325 vj[1]: 0.211325 x[0]: 0.211325 y[0]: 1.21132 JxW[0]: 0.250000 distance = 1.00000 +DEAL::vi[0]: 0.788675 vj[0]: 0.211325 x[1]: 0.211325 y[1]: 1.78868 JxW[1]: 0.250000 distance = 1.57735 +DEAL::vi[0]: 0.788675 vj[1]: 0.788675 x[1]: 0.211325 y[1]: 1.78868 JxW[1]: 0.250000 distance = 1.57735 +DEAL::vi[1]: 0.211325 vj[0]: 0.211325 x[1]: 0.211325 y[1]: 1.78868 JxW[1]: 0.250000 distance = 1.57735 +DEAL::vi[1]: 0.211325 vj[1]: 0.788675 x[1]: 0.211325 y[1]: 1.78868 JxW[1]: 0.250000 distance = 1.57735 +DEAL::vi[0]: 0.211325 vj[0]: 0.788675 x[2]: 0.788675 y[2]: 1.21132 JxW[2]: 0.250000 distance = 0.422650 +DEAL::vi[0]: 0.211325 vj[1]: 0.211325 x[2]: 0.788675 y[2]: 1.21132 JxW[2]: 0.250000 distance = 0.422650 +DEAL::vi[1]: 0.788675 vj[0]: 0.788675 x[2]: 0.788675 y[2]: 1.21132 JxW[2]: 0.250000 distance = 0.422650 +DEAL::vi[1]: 0.788675 vj[1]: 0.211325 x[2]: 0.788675 y[2]: 1.21132 JxW[2]: 0.250000 distance = 0.422650 +DEAL::vi[0]: 0.211325 vj[0]: 0.211325 x[3]: 0.788675 y[3]: 1.78868 JxW[3]: 0.250000 distance = 1.00000 +DEAL::vi[0]: 0.211325 vj[1]: 0.788675 x[3]: 0.788675 y[3]: 1.78868 JxW[3]: 0.250000 distance = 1.00000 +DEAL::vi[1]: 0.788675 vj[0]: 0.211325 x[3]: 0.788675 y[3]: 1.78868 JxW[3]: 0.250000 distance = 1.00000 +DEAL::vi[1]: 0.788675 vj[1]: 0.788675 x[3]: 0.788675 y[3]: 1.78868 JxW[3]: 0.250000 distance = 1.00000 +DEAL::matrix = +DEAL::0.250000 0.333333 +DEAL::0.166667 0.250000 diff --git a/tests/fecoupling/fe_coupling_03.cc b/tests/fecoupling/fe_coupling_03.cc new file mode 100644 index 0000000000..bdc4f2bc8d --- /dev/null +++ b/tests/fecoupling/fe_coupling_03.cc @@ -0,0 +1,169 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test QuadratureCouplingType::matching integration in FECouplingValues +// with DoFCouplingType::contiguous. Typical for DG interfaces. + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + +#include "../test_grids.h" + +int +main() +{ + initlog(); + constexpr unsigned int dim = 2; + + Triangulation tria; + TestGrids::hyper_line(tria, 2); + + DoFHandler dofh(tria); + FE_DGQ fe(1); + dofh.distribute_dofs(fe); + + UpdateFlags update_flags = + update_quadrature_points | update_values | update_JxW_values; + + QGauss quad(fe.degree + 1); + + FEFaceValues fv1(fe, quad, update_flags); + FEFaceValues fv2(fe, quad, update_flags); + FEInterfaceValues fiv(fe, quad, update_flags); + + std::vector v1(fe.dofs_per_cell); + std::vector 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 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 matrix(fcv.n_coupling_dofs(), fcv.n_coupling_dofs()); + FEValuesExtractors::Scalar scalar(0); + + const auto first = fcv.get_first_extractor(scalar); + const auto second = fcv.get_second_extractor(scalar); + + // We need to loop over all the coupling quadrature points + for (const auto q : fcv.quadrature_point_indices()) + { + const auto &[x, y] = fcv.quadrature_point(q); + for (const auto i : fcv.coupling_dof_indices()) + { + const auto &lvi = fcv[first].value(i, q); + const auto &rvi = fcv[second].value(i, q); + const auto jump_i = lvi - rvi; + const auto &jump_ifiv = fiv[scalar].jump_in_values(i, q); + + + for (const auto j : fcv.coupling_dof_indices()) + { + const auto &lvj = fcv[first].value(j, q); + const auto &rvj = fcv[second].value(j, q); + + const auto jump_j = lvj - rvj; + const auto &jump_jfiv = fiv[scalar].jump_in_values(j, q); + + deallog << std::left << "fcv -- " // + << "[vi][" << i << "]: " << std::setw(10) << jump_i // + << "[vj][" << j << "]: " << std::setw(10) << jump_j // + << "x[" << q << "]: " << std::setw(10) << x // + << " y[" << q << "]: " << std::setw(10) << y // + << "JxW[" << q << "]: " << std::setw(10) << fcv.JxW(q) + << std::endl; + deallog << std::left << "fiv -- " // + << "[vi][" << i << "]: " << std::setw(10) << jump_ifiv // + << "[vj][" << j << "]: " << std::setw(10) << jump_jfiv // + << std::endl; + matrix(i, j) += jump_i * jump_j * fcv.JxW(q); + } + } + } + // Now print the matrix: + deallog << "matrix = " << std::endl; + for (unsigned int i = 0; i < matrix.m(); ++i) + { + for (unsigned int j = 0; j < matrix.n(); ++j) + { + deallog << matrix(i, j) << " "; + } + deallog << std::endl; + } +} diff --git a/tests/fecoupling/fe_coupling_03.output b/tests/fecoupling/fe_coupling_03.output new file mode 100644 index 0000000000..14ce313dd2 --- /dev/null +++ b/tests/fecoupling/fe_coupling_03.output @@ -0,0 +1,271 @@ + +DEAL::dofs_per_cell = 4 +DEAL::n_q_points_per_cell = 2 +DEAL::n_q_points = 2 +DEAL::n_interface_dofs = 8 +DEAL::n_coupling_dofs = 8 +DEAL::fcv -- [vi][0]: -0.00000 [vj][0]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][0]: 0.00000 +DEAL::fcv -- [vi][0]: -0.00000 [vj][1]: 0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][1]: 0.788675 +DEAL::fcv -- [vi][0]: -0.00000 [vj][2]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][2]: 0.00000 +DEAL::fcv -- [vi][0]: -0.00000 [vj][3]: 0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][3]: 0.211325 +DEAL::fcv -- [vi][0]: -0.00000 [vj][4]: -0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][4]: -0.788675 +DEAL::fcv -- [vi][0]: -0.00000 [vj][5]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][5]: 0.00000 +DEAL::fcv -- [vi][0]: -0.00000 [vj][6]: -0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][6]: -0.211325 +DEAL::fcv -- [vi][0]: -0.00000 [vj][7]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][7]: 0.00000 +DEAL::fcv -- [vi][1]: 0.788675 [vj][0]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][1]: 0.788675 [vj][0]: 0.00000 +DEAL::fcv -- [vi][1]: 0.788675 [vj][1]: 0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][1]: 0.788675 [vj][1]: 0.788675 +DEAL::fcv -- [vi][1]: 0.788675 [vj][2]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][1]: 0.788675 [vj][2]: 0.00000 +DEAL::fcv -- [vi][1]: 0.788675 [vj][3]: 0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][1]: 0.788675 [vj][3]: 0.211325 +DEAL::fcv -- [vi][1]: 0.788675 [vj][4]: -0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][1]: 0.788675 [vj][4]: -0.788675 +DEAL::fcv -- [vi][1]: 0.788675 [vj][5]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][1]: 0.788675 [vj][5]: 0.00000 +DEAL::fcv -- [vi][1]: 0.788675 [vj][6]: -0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][1]: 0.788675 [vj][6]: -0.211325 +DEAL::fcv -- [vi][1]: 0.788675 [vj][7]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][1]: 0.788675 [vj][7]: 0.00000 +DEAL::fcv -- [vi][2]: -0.00000 [vj][0]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][0]: 0.00000 +DEAL::fcv -- [vi][2]: -0.00000 [vj][1]: 0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][1]: 0.788675 +DEAL::fcv -- [vi][2]: -0.00000 [vj][2]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][2]: 0.00000 +DEAL::fcv -- [vi][2]: -0.00000 [vj][3]: 0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][3]: 0.211325 +DEAL::fcv -- [vi][2]: -0.00000 [vj][4]: -0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][4]: -0.788675 +DEAL::fcv -- [vi][2]: -0.00000 [vj][5]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][5]: 0.00000 +DEAL::fcv -- [vi][2]: -0.00000 [vj][6]: -0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][6]: -0.211325 +DEAL::fcv -- [vi][2]: -0.00000 [vj][7]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][7]: 0.00000 +DEAL::fcv -- [vi][3]: 0.211325 [vj][0]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][3]: 0.211325 [vj][0]: 0.00000 +DEAL::fcv -- [vi][3]: 0.211325 [vj][1]: 0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][3]: 0.211325 [vj][1]: 0.788675 +DEAL::fcv -- [vi][3]: 0.211325 [vj][2]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][3]: 0.211325 [vj][2]: 0.00000 +DEAL::fcv -- [vi][3]: 0.211325 [vj][3]: 0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][3]: 0.211325 [vj][3]: 0.211325 +DEAL::fcv -- [vi][3]: 0.211325 [vj][4]: -0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][3]: 0.211325 [vj][4]: -0.788675 +DEAL::fcv -- [vi][3]: 0.211325 [vj][5]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][3]: 0.211325 [vj][5]: 0.00000 +DEAL::fcv -- [vi][3]: 0.211325 [vj][6]: -0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][3]: 0.211325 [vj][6]: -0.211325 +DEAL::fcv -- [vi][3]: 0.211325 [vj][7]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][3]: 0.211325 [vj][7]: 0.00000 +DEAL::fcv -- [vi][4]: -0.788675 [vj][0]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][4]: -0.788675 [vj][0]: 0.00000 +DEAL::fcv -- [vi][4]: -0.788675 [vj][1]: 0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][4]: -0.788675 [vj][1]: 0.788675 +DEAL::fcv -- [vi][4]: -0.788675 [vj][2]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][4]: -0.788675 [vj][2]: 0.00000 +DEAL::fcv -- [vi][4]: -0.788675 [vj][3]: 0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][4]: -0.788675 [vj][3]: 0.211325 +DEAL::fcv -- [vi][4]: -0.788675 [vj][4]: -0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][4]: -0.788675 [vj][4]: -0.788675 +DEAL::fcv -- [vi][4]: -0.788675 [vj][5]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][4]: -0.788675 [vj][5]: 0.00000 +DEAL::fcv -- [vi][4]: -0.788675 [vj][6]: -0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][4]: -0.788675 [vj][6]: -0.211325 +DEAL::fcv -- [vi][4]: -0.788675 [vj][7]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][4]: -0.788675 [vj][7]: 0.00000 +DEAL::fcv -- [vi][5]: 0.00000 [vj][0]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][0]: 0.00000 +DEAL::fcv -- [vi][5]: 0.00000 [vj][1]: 0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][1]: 0.788675 +DEAL::fcv -- [vi][5]: 0.00000 [vj][2]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][2]: 0.00000 +DEAL::fcv -- [vi][5]: 0.00000 [vj][3]: 0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][3]: 0.211325 +DEAL::fcv -- [vi][5]: 0.00000 [vj][4]: -0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][4]: -0.788675 +DEAL::fcv -- [vi][5]: 0.00000 [vj][5]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][5]: 0.00000 +DEAL::fcv -- [vi][5]: 0.00000 [vj][6]: -0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][6]: -0.211325 +DEAL::fcv -- [vi][5]: 0.00000 [vj][7]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][7]: 0.00000 +DEAL::fcv -- [vi][6]: -0.211325 [vj][0]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][6]: -0.211325 [vj][0]: 0.00000 +DEAL::fcv -- [vi][6]: -0.211325 [vj][1]: 0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][6]: -0.211325 [vj][1]: 0.788675 +DEAL::fcv -- [vi][6]: -0.211325 [vj][2]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][6]: -0.211325 [vj][2]: 0.00000 +DEAL::fcv -- [vi][6]: -0.211325 [vj][3]: 0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][6]: -0.211325 [vj][3]: 0.211325 +DEAL::fcv -- [vi][6]: -0.211325 [vj][4]: -0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][6]: -0.211325 [vj][4]: -0.788675 +DEAL::fcv -- [vi][6]: -0.211325 [vj][5]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][6]: -0.211325 [vj][5]: 0.00000 +DEAL::fcv -- [vi][6]: -0.211325 [vj][6]: -0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][6]: -0.211325 [vj][6]: -0.211325 +DEAL::fcv -- [vi][6]: -0.211325 [vj][7]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][6]: -0.211325 [vj][7]: 0.00000 +DEAL::fcv -- [vi][7]: 0.00000 [vj][0]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][0]: 0.00000 +DEAL::fcv -- [vi][7]: 0.00000 [vj][1]: 0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][1]: 0.788675 +DEAL::fcv -- [vi][7]: 0.00000 [vj][2]: -0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][2]: 0.00000 +DEAL::fcv -- [vi][7]: 0.00000 [vj][3]: 0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][3]: 0.211325 +DEAL::fcv -- [vi][7]: 0.00000 [vj][4]: -0.788675 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][4]: -0.788675 +DEAL::fcv -- [vi][7]: 0.00000 [vj][5]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][5]: 0.00000 +DEAL::fcv -- [vi][7]: 0.00000 [vj][6]: -0.211325 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][6]: -0.211325 +DEAL::fcv -- [vi][7]: 0.00000 [vj][7]: 0.00000 x[0]: 1.00000 0.211325 y[0]: 1.00000 0.211325JxW[0]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][7]: 0.00000 +DEAL::fcv -- [vi][0]: -0.00000 [vj][0]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][0]: 0.00000 +DEAL::fcv -- [vi][0]: -0.00000 [vj][1]: 0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][1]: 0.211325 +DEAL::fcv -- [vi][0]: -0.00000 [vj][2]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][2]: 0.00000 +DEAL::fcv -- [vi][0]: -0.00000 [vj][3]: 0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][3]: 0.788675 +DEAL::fcv -- [vi][0]: -0.00000 [vj][4]: -0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][4]: -0.211325 +DEAL::fcv -- [vi][0]: -0.00000 [vj][5]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][5]: 0.00000 +DEAL::fcv -- [vi][0]: -0.00000 [vj][6]: -0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][6]: -0.788675 +DEAL::fcv -- [vi][0]: -0.00000 [vj][7]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][0]: 0.00000 [vj][7]: 0.00000 +DEAL::fcv -- [vi][1]: 0.211325 [vj][0]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][1]: 0.211325 [vj][0]: 0.00000 +DEAL::fcv -- [vi][1]: 0.211325 [vj][1]: 0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][1]: 0.211325 [vj][1]: 0.211325 +DEAL::fcv -- [vi][1]: 0.211325 [vj][2]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][1]: 0.211325 [vj][2]: 0.00000 +DEAL::fcv -- [vi][1]: 0.211325 [vj][3]: 0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][1]: 0.211325 [vj][3]: 0.788675 +DEAL::fcv -- [vi][1]: 0.211325 [vj][4]: -0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][1]: 0.211325 [vj][4]: -0.211325 +DEAL::fcv -- [vi][1]: 0.211325 [vj][5]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][1]: 0.211325 [vj][5]: 0.00000 +DEAL::fcv -- [vi][1]: 0.211325 [vj][6]: -0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][1]: 0.211325 [vj][6]: -0.788675 +DEAL::fcv -- [vi][1]: 0.211325 [vj][7]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][1]: 0.211325 [vj][7]: 0.00000 +DEAL::fcv -- [vi][2]: -0.00000 [vj][0]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][0]: 0.00000 +DEAL::fcv -- [vi][2]: -0.00000 [vj][1]: 0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][1]: 0.211325 +DEAL::fcv -- [vi][2]: -0.00000 [vj][2]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][2]: 0.00000 +DEAL::fcv -- [vi][2]: -0.00000 [vj][3]: 0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][3]: 0.788675 +DEAL::fcv -- [vi][2]: -0.00000 [vj][4]: -0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][4]: -0.211325 +DEAL::fcv -- [vi][2]: -0.00000 [vj][5]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][5]: 0.00000 +DEAL::fcv -- [vi][2]: -0.00000 [vj][6]: -0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][6]: -0.788675 +DEAL::fcv -- [vi][2]: -0.00000 [vj][7]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][2]: 0.00000 [vj][7]: 0.00000 +DEAL::fcv -- [vi][3]: 0.788675 [vj][0]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][3]: 0.788675 [vj][0]: 0.00000 +DEAL::fcv -- [vi][3]: 0.788675 [vj][1]: 0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][3]: 0.788675 [vj][1]: 0.211325 +DEAL::fcv -- [vi][3]: 0.788675 [vj][2]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][3]: 0.788675 [vj][2]: 0.00000 +DEAL::fcv -- [vi][3]: 0.788675 [vj][3]: 0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][3]: 0.788675 [vj][3]: 0.788675 +DEAL::fcv -- [vi][3]: 0.788675 [vj][4]: -0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][3]: 0.788675 [vj][4]: -0.211325 +DEAL::fcv -- [vi][3]: 0.788675 [vj][5]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][3]: 0.788675 [vj][5]: 0.00000 +DEAL::fcv -- [vi][3]: 0.788675 [vj][6]: -0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][3]: 0.788675 [vj][6]: -0.788675 +DEAL::fcv -- [vi][3]: 0.788675 [vj][7]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][3]: 0.788675 [vj][7]: 0.00000 +DEAL::fcv -- [vi][4]: -0.211325 [vj][0]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][4]: -0.211325 [vj][0]: 0.00000 +DEAL::fcv -- [vi][4]: -0.211325 [vj][1]: 0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][4]: -0.211325 [vj][1]: 0.211325 +DEAL::fcv -- [vi][4]: -0.211325 [vj][2]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][4]: -0.211325 [vj][2]: 0.00000 +DEAL::fcv -- [vi][4]: -0.211325 [vj][3]: 0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][4]: -0.211325 [vj][3]: 0.788675 +DEAL::fcv -- [vi][4]: -0.211325 [vj][4]: -0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][4]: -0.211325 [vj][4]: -0.211325 +DEAL::fcv -- [vi][4]: -0.211325 [vj][5]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][4]: -0.211325 [vj][5]: 0.00000 +DEAL::fcv -- [vi][4]: -0.211325 [vj][6]: -0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][4]: -0.211325 [vj][6]: -0.788675 +DEAL::fcv -- [vi][4]: -0.211325 [vj][7]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][4]: -0.211325 [vj][7]: 0.00000 +DEAL::fcv -- [vi][5]: 0.00000 [vj][0]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][0]: 0.00000 +DEAL::fcv -- [vi][5]: 0.00000 [vj][1]: 0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][1]: 0.211325 +DEAL::fcv -- [vi][5]: 0.00000 [vj][2]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][2]: 0.00000 +DEAL::fcv -- [vi][5]: 0.00000 [vj][3]: 0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][3]: 0.788675 +DEAL::fcv -- [vi][5]: 0.00000 [vj][4]: -0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][4]: -0.211325 +DEAL::fcv -- [vi][5]: 0.00000 [vj][5]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][5]: 0.00000 +DEAL::fcv -- [vi][5]: 0.00000 [vj][6]: -0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][6]: -0.788675 +DEAL::fcv -- [vi][5]: 0.00000 [vj][7]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][5]: 0.00000 [vj][7]: 0.00000 +DEAL::fcv -- [vi][6]: -0.788675 [vj][0]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][6]: -0.788675 [vj][0]: 0.00000 +DEAL::fcv -- [vi][6]: -0.788675 [vj][1]: 0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][6]: -0.788675 [vj][1]: 0.211325 +DEAL::fcv -- [vi][6]: -0.788675 [vj][2]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][6]: -0.788675 [vj][2]: 0.00000 +DEAL::fcv -- [vi][6]: -0.788675 [vj][3]: 0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][6]: -0.788675 [vj][3]: 0.788675 +DEAL::fcv -- [vi][6]: -0.788675 [vj][4]: -0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][6]: -0.788675 [vj][4]: -0.211325 +DEAL::fcv -- [vi][6]: -0.788675 [vj][5]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][6]: -0.788675 [vj][5]: 0.00000 +DEAL::fcv -- [vi][6]: -0.788675 [vj][6]: -0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][6]: -0.788675 [vj][6]: -0.788675 +DEAL::fcv -- [vi][6]: -0.788675 [vj][7]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][6]: -0.788675 [vj][7]: 0.00000 +DEAL::fcv -- [vi][7]: 0.00000 [vj][0]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][0]: 0.00000 +DEAL::fcv -- [vi][7]: 0.00000 [vj][1]: 0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][1]: 0.211325 +DEAL::fcv -- [vi][7]: 0.00000 [vj][2]: -0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][2]: 0.00000 +DEAL::fcv -- [vi][7]: 0.00000 [vj][3]: 0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][3]: 0.788675 +DEAL::fcv -- [vi][7]: 0.00000 [vj][4]: -0.211325 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][4]: -0.211325 +DEAL::fcv -- [vi][7]: 0.00000 [vj][5]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][5]: 0.00000 +DEAL::fcv -- [vi][7]: 0.00000 [vj][6]: -0.788675 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][6]: -0.788675 +DEAL::fcv -- [vi][7]: 0.00000 [vj][7]: 0.00000 x[1]: 1.00000 0.788675 y[1]: 1.00000 0.788675JxW[1]: 0.500000 +DEAL::fiv -- [vi][7]: 0.00000 [vj][7]: 0.00000 +DEAL::matrix = +DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL::0.00000 0.333333 0.00000 0.166667 -0.333333 0.00000 -0.166667 0.00000 +DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL::0.00000 0.166667 0.00000 0.333333 -0.166667 0.00000 -0.333333 0.00000 +DEAL::0.00000 -0.333333 0.00000 -0.166667 0.333333 0.00000 0.166667 0.00000 +DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL::0.00000 -0.166667 0.00000 -0.333333 0.166667 0.00000 0.333333 0.00000 +DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 diff --git a/tests/fecoupling/fe_coupling_04.cc b/tests/fecoupling/fe_coupling_04.cc new file mode 100644 index 0000000000..2342960c9c --- /dev/null +++ b/tests/fecoupling/fe_coupling_04.cc @@ -0,0 +1,154 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test QuadratureCouplingType::reorder integration in FECouplingValues +// with DoFCouplingType::contiguous. Typical for bulk-surface coupling. + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + +#include "../test_grids.h" + +int +main() +{ + initlog(); + constexpr unsigned int dim = 2; + constexpr unsigned int codim = dim - 1; + + Triangulation tria; + TestGrids::hyper_line(tria, 2); + + Triangulation codim_tria; + TestGrids::hyper_line(codim_tria, 2); + + DoFHandler dofh(tria); + FE_Q fe(1); + dofh.distribute_dofs(fe); + + DoFHandler codim_dofh(codim_tria); + FE_Q codim_fe(1); + codim_dofh.distribute_dofs(codim_fe); + + UpdateFlags update_flags = + update_quadrature_points | update_values | update_JxW_values; + + QGauss quad(fe.degree + 1); + + FEFaceValues fv1(fe, quad, update_flags); + FEValues fv2(codim_fe, quad, update_flags); + + std::vector v1(fe.dofs_per_cell); + std::vector 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 fcv(fv1, + fv2, + DoFCouplingType::contiguous, + QuadratureCouplingType::reorder); + + FullMatrix matrix(fcv.n_coupling_dofs(), fcv.n_coupling_dofs()); + FEValuesExtractors::Scalar scalar(0); + + const auto bulk = fcv.get_first_extractor(scalar); + const auto surface = fcv.get_second_extractor(scalar); + + deallog << "Renumbering vectors: " << std::endl; + + for (const auto i : fcv.coupling_dof_indices()) + { + auto id_first = fcv.coupling_dof_to_dof_indices(i).first; + auto id_second = fcv.coupling_dof_to_dof_indices(i).second; + deallog << "i: " << i << " bulk: " << (int)id_first + << " surface: " << (int)id_second << std::endl; + } + + // We need to loop over all the coupling quadrature points + for (const auto q : fcv.quadrature_point_indices()) + { + const auto &[x, y] = fcv.quadrature_point(q); + for (const auto i : fcv.coupling_dof_indices()) + { + const auto &bulk_vi = fcv[bulk].value(i, q); + const auto &surface_vi = fcv[surface].value(i, q); + + for (const auto j : fcv.coupling_dof_indices()) + { + const auto &bulk_vj = fcv[bulk].value(j, q); + const auto &surface_vj = fcv[surface].value(j, q); + + deallog << std::left // + << "bvi[" << i << "]: " << std::setw(10) << bulk_vi // + << "bvj[" << j << "]: " << std::setw(10) << bulk_vj // + << "svi[" << i << "]: " << std::setw(10) << surface_vi // + << "svj[" << j << "]: " << std::setw(10) << surface_vj // + << "x[" << q << "]: " << std::setw(10) << x // + << " y[" << q << "]: " << std::setw(10) << y // + << std::endl; + matrix(i, j) += (bulk_vi * bulk_vj + // + 10 * bulk_vi * surface_vj + // + 100 * surface_vi * bulk_vj + // + 1000 * surface_vi * surface_vj) // + * fcv.JxW(q); + } + } + } + // Now print the matrix: + deallog << "matrix = " << std::endl << std::left; + for (unsigned int i = 0; i < matrix.m(); ++i) + { + for (unsigned int j = 0; j < matrix.n(); ++j) + { + deallog << std::setw(8) << matrix(i, j) << " "; + } + deallog << std::endl; + } +} diff --git a/tests/fecoupling/fe_coupling_04.output b/tests/fecoupling/fe_coupling_04.output new file mode 100644 index 0000000000..b182b41046 --- /dev/null +++ b/tests/fecoupling/fe_coupling_04.output @@ -0,0 +1,88 @@ + +DEAL::Bulk dofs: 4, Surface dofs: 2 +DEAL::Renumbering vectors: +DEAL::i: 0 bulk: 0 surface: -1 +DEAL::i: 1 bulk: 1 surface: -1 +DEAL::i: 2 bulk: 2 surface: -1 +DEAL::i: 3 bulk: 3 surface: -1 +DEAL::i: 4 bulk: -1 surface: 0 +DEAL::i: 5 bulk: -1 surface: 1 +DEAL::bvi[0]: 0.788675 bvj[0]: 0.788675 svi[0]: 0.00000 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.788675 bvj[1]: 0.211325 svi[0]: 0.00000 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.788675 bvj[2]: 0.00000 svi[0]: 0.00000 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.788675 bvj[3]: 0.00000 svi[0]: 0.00000 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.788675 bvj[4]: 0.00000 svi[0]: 0.00000 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.788675 bvj[5]: 0.00000 svi[0]: 0.00000 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[0]: 0.788675 svi[1]: 0.00000 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[1]: 0.211325 svi[1]: 0.00000 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[2]: 0.00000 svi[1]: 0.00000 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[3]: 0.00000 svi[1]: 0.00000 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[4]: 0.00000 svi[1]: 0.00000 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[5]: 0.00000 svi[1]: 0.00000 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[0]: 0.788675 svi[2]: 0.00000 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[1]: 0.211325 svi[2]: 0.00000 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[2]: 0.00000 svi[2]: 0.00000 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[3]: 0.00000 svi[2]: 0.00000 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[4]: 0.00000 svi[2]: 0.00000 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[5]: 0.00000 svi[2]: 0.00000 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[0]: 0.788675 svi[3]: 0.00000 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[1]: 0.211325 svi[3]: 0.00000 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[2]: 0.00000 svi[3]: 0.00000 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[3]: 0.00000 svi[3]: 0.00000 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[4]: 0.00000 svi[3]: 0.00000 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[5]: 0.00000 svi[3]: 0.00000 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[0]: 0.788675 svi[4]: 0.788675 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[1]: 0.211325 svi[4]: 0.788675 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[2]: 0.00000 svi[4]: 0.788675 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[3]: 0.00000 svi[4]: 0.788675 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[4]: 0.00000 svi[4]: 0.788675 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[5]: 0.00000 svi[4]: 0.788675 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[0]: 0.788675 svi[5]: 0.211325 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[1]: 0.211325 svi[5]: 0.211325 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[2]: 0.00000 svi[5]: 0.211325 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[3]: 0.00000 svi[5]: 0.211325 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[4]: 0.00000 svi[5]: 0.211325 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[5]: 0.00000 svi[5]: 0.211325 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.211325 bvj[0]: 0.211325 svi[0]: 0.00000 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[0]: 0.211325 bvj[1]: 0.788675 svi[0]: 0.00000 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[0]: 0.211325 bvj[2]: 0.00000 svi[0]: 0.00000 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[0]: 0.211325 bvj[3]: 0.00000 svi[0]: 0.00000 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[0]: 0.211325 bvj[4]: 0.00000 svi[0]: 0.00000 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[0]: 0.211325 bvj[5]: 0.00000 svi[0]: 0.00000 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[0]: 0.211325 svi[1]: 0.00000 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[1]: 0.788675 svi[1]: 0.00000 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[2]: 0.00000 svi[1]: 0.00000 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[3]: 0.00000 svi[1]: 0.00000 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[4]: 0.00000 svi[1]: 0.00000 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[5]: 0.00000 svi[1]: 0.00000 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[0]: 0.211325 svi[2]: 0.00000 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[1]: 0.788675 svi[2]: 0.00000 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[2]: 0.00000 svi[2]: 0.00000 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[3]: 0.00000 svi[2]: 0.00000 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[4]: 0.00000 svi[2]: 0.00000 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[5]: 0.00000 svi[2]: 0.00000 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[0]: 0.211325 svi[3]: 0.00000 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[1]: 0.788675 svi[3]: 0.00000 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[2]: 0.00000 svi[3]: 0.00000 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[3]: 0.00000 svi[3]: 0.00000 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[4]: 0.00000 svi[3]: 0.00000 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[5]: 0.00000 svi[3]: 0.00000 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[0]: 0.211325 svi[4]: 0.211325 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[1]: 0.788675 svi[4]: 0.211325 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[2]: 0.00000 svi[4]: 0.211325 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[3]: 0.00000 svi[4]: 0.211325 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[4]: 0.00000 svi[4]: 0.211325 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[5]: 0.00000 svi[4]: 0.211325 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[0]: 0.211325 svi[5]: 0.788675 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[1]: 0.788675 svi[5]: 0.788675 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[2]: 0.00000 svi[5]: 0.788675 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[3]: 0.00000 svi[5]: 0.788675 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[4]: 0.00000 svi[5]: 0.788675 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[5]: 0.00000 svi[5]: 0.788675 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::matrix = +DEAL::0.333333 0.166667 0.00000 0.00000 3.33333 1.66667 +DEAL::0.166667 0.333333 0.00000 0.00000 1.66667 3.33333 +DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL::33.3333 16.6667 0.00000 0.00000 333.333 166.667 +DEAL::16.6667 33.3333 0.00000 0.00000 166.667 333.333 diff --git a/tests/fecoupling/fe_coupling_05.cc b/tests/fecoupling/fe_coupling_05.cc new file mode 100644 index 0000000000..b83e4289ef --- /dev/null +++ b/tests/fecoupling/fe_coupling_05.cc @@ -0,0 +1,154 @@ +// --------------------------------------------------------------------- +// +// Copyright (C) 2023 by the deal.II authors +// +// This file is part of the deal.II library. +// +// The deal.II library is free software; you can use it, redistribute +// it, and/or modify it under the terms of the GNU Lesser General +// Public License as published by the Free Software Foundation; either +// version 2.1 of the License, or (at your option) any later version. +// The full text of the license can be found in the file LICENSE.md at +// the top level directory of deal.II. +// +// --------------------------------------------------------------------- + + +// Test QuadratureCouplingType::matching integration in FECouplingValues +// with DoFCouplingType::contiguous. Typical for bulk-surface coupling. + +#include + +#include +#include +#include +#include +#include +#include + +#include +#include +#include + +#include +#include +#include + +#include "../tests.h" + +#include "../test_grids.h" + +int +main() +{ + initlog(); + constexpr unsigned int dim = 2; + constexpr unsigned int codim = dim - 1; + + Triangulation tria; + TestGrids::hyper_line(tria, 2); + + Triangulation codim_tria; + TestGrids::hyper_line(codim_tria, 2); + + DoFHandler dofh(tria); + FE_Q fe(1); + dofh.distribute_dofs(fe); + + DoFHandler codim_dofh(codim_tria); + FE_Q codim_fe(1); + codim_dofh.distribute_dofs(codim_fe); + + UpdateFlags update_flags = + update_quadrature_points | update_values | update_JxW_values; + + QGauss quad(fe.degree + 1); + + FEFaceValues fv1(fe, quad, update_flags); + FEValues fv2(codim_fe, quad, update_flags); + + std::vector v1(fe.dofs_per_cell); + std::vector 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 fcv(fv1, + fv2, + DoFCouplingType::contiguous, + QuadratureCouplingType::matching); + + FullMatrix matrix(fcv.n_coupling_dofs(), fcv.n_coupling_dofs()); + FEValuesExtractors::Scalar scalar(0); + + const auto bulk = fcv.get_first_extractor(scalar); + const auto surface = fcv.get_second_extractor(scalar); + + deallog << "Renumbering vectors: " << std::endl; + + for (const auto i : fcv.coupling_dof_indices()) + { + auto id_first = fcv.coupling_dof_to_dof_indices(i).first; + auto id_second = fcv.coupling_dof_to_dof_indices(i).second; + deallog << "i: " << i << " bulk: " << (int)id_first + << " surface: " << (int)id_second << std::endl; + } + + // We need to loop over all the coupling quadrature points + for (const auto q : fcv.quadrature_point_indices()) + { + const auto &[x, y] = fcv.quadrature_point(q); + for (const auto i : fcv.coupling_dof_indices()) + { + const auto &bulk_vi = fcv[bulk].value(i, q); + const auto &surface_vi = fcv[surface].value(i, q); + + for (const auto j : fcv.coupling_dof_indices()) + { + const auto &bulk_vj = fcv[bulk].value(j, q); + const auto &surface_vj = fcv[surface].value(j, q); + + deallog << std::left // + << "bvi[" << i << "]: " << std::setw(10) << bulk_vi // + << "bvj[" << j << "]: " << std::setw(10) << bulk_vj // + << "svi[" << i << "]: " << std::setw(10) << surface_vi // + << "svj[" << j << "]: " << std::setw(10) << surface_vj // + << "x[" << q << "]: " << std::setw(10) << x // + << " y[" << q << "]: " << std::setw(10) << y // + << std::endl; + matrix(i, j) += (bulk_vi * bulk_vj + // + 10 * bulk_vi * surface_vj + // + 100 * surface_vi * bulk_vj + // + 1000 * surface_vi * surface_vj) // + * fcv.JxW(q); + } + } + } + // Now print the matrix: + deallog << "matrix = " << std::endl << std::left; + for (unsigned int i = 0; i < matrix.m(); ++i) + { + for (unsigned int j = 0; j < matrix.n(); ++j) + { + deallog << std::setw(8) << matrix(i, j) << " "; + } + deallog << std::endl; + } +} diff --git a/tests/fecoupling/fe_coupling_05.output b/tests/fecoupling/fe_coupling_05.output new file mode 100644 index 0000000000..b182b41046 --- /dev/null +++ b/tests/fecoupling/fe_coupling_05.output @@ -0,0 +1,88 @@ + +DEAL::Bulk dofs: 4, Surface dofs: 2 +DEAL::Renumbering vectors: +DEAL::i: 0 bulk: 0 surface: -1 +DEAL::i: 1 bulk: 1 surface: -1 +DEAL::i: 2 bulk: 2 surface: -1 +DEAL::i: 3 bulk: 3 surface: -1 +DEAL::i: 4 bulk: -1 surface: 0 +DEAL::i: 5 bulk: -1 surface: 1 +DEAL::bvi[0]: 0.788675 bvj[0]: 0.788675 svi[0]: 0.00000 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.788675 bvj[1]: 0.211325 svi[0]: 0.00000 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.788675 bvj[2]: 0.00000 svi[0]: 0.00000 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.788675 bvj[3]: 0.00000 svi[0]: 0.00000 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.788675 bvj[4]: 0.00000 svi[0]: 0.00000 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.788675 bvj[5]: 0.00000 svi[0]: 0.00000 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[0]: 0.788675 svi[1]: 0.00000 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[1]: 0.211325 svi[1]: 0.00000 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[2]: 0.00000 svi[1]: 0.00000 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[3]: 0.00000 svi[1]: 0.00000 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[4]: 0.00000 svi[1]: 0.00000 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[1]: 0.211325 bvj[5]: 0.00000 svi[1]: 0.00000 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[0]: 0.788675 svi[2]: 0.00000 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[1]: 0.211325 svi[2]: 0.00000 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[2]: 0.00000 svi[2]: 0.00000 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[3]: 0.00000 svi[2]: 0.00000 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[4]: 0.00000 svi[2]: 0.00000 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[2]: 0.00000 bvj[5]: 0.00000 svi[2]: 0.00000 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[0]: 0.788675 svi[3]: 0.00000 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[1]: 0.211325 svi[3]: 0.00000 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[2]: 0.00000 svi[3]: 0.00000 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[3]: 0.00000 svi[3]: 0.00000 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[4]: 0.00000 svi[3]: 0.00000 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[3]: 0.00000 bvj[5]: 0.00000 svi[3]: 0.00000 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[0]: 0.788675 svi[4]: 0.788675 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[1]: 0.211325 svi[4]: 0.788675 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[2]: 0.00000 svi[4]: 0.788675 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[3]: 0.00000 svi[4]: 0.788675 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[4]: 0.00000 svi[4]: 0.788675 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[4]: 0.00000 bvj[5]: 0.00000 svi[4]: 0.788675 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[0]: 0.788675 svi[5]: 0.211325 svj[0]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[1]: 0.211325 svi[5]: 0.211325 svj[1]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[2]: 0.00000 svi[5]: 0.211325 svj[2]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[3]: 0.00000 svi[5]: 0.211325 svj[3]: 0.00000 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[4]: 0.00000 svi[5]: 0.211325 svj[4]: 0.788675 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[5]: 0.00000 bvj[5]: 0.00000 svi[5]: 0.211325 svj[5]: 0.211325 x[0]: 0.211325 0.00000 y[0]: 0.211325 0.00000 +DEAL::bvi[0]: 0.211325 bvj[0]: 0.211325 svi[0]: 0.00000 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[0]: 0.211325 bvj[1]: 0.788675 svi[0]: 0.00000 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[0]: 0.211325 bvj[2]: 0.00000 svi[0]: 0.00000 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[0]: 0.211325 bvj[3]: 0.00000 svi[0]: 0.00000 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[0]: 0.211325 bvj[4]: 0.00000 svi[0]: 0.00000 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[0]: 0.211325 bvj[5]: 0.00000 svi[0]: 0.00000 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[0]: 0.211325 svi[1]: 0.00000 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[1]: 0.788675 svi[1]: 0.00000 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[2]: 0.00000 svi[1]: 0.00000 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[3]: 0.00000 svi[1]: 0.00000 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[4]: 0.00000 svi[1]: 0.00000 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[1]: 0.788675 bvj[5]: 0.00000 svi[1]: 0.00000 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[0]: 0.211325 svi[2]: 0.00000 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[1]: 0.788675 svi[2]: 0.00000 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[2]: 0.00000 svi[2]: 0.00000 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[3]: 0.00000 svi[2]: 0.00000 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[4]: 0.00000 svi[2]: 0.00000 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[2]: 0.00000 bvj[5]: 0.00000 svi[2]: 0.00000 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[0]: 0.211325 svi[3]: 0.00000 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[1]: 0.788675 svi[3]: 0.00000 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[2]: 0.00000 svi[3]: 0.00000 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[3]: 0.00000 svi[3]: 0.00000 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[4]: 0.00000 svi[3]: 0.00000 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[3]: 0.00000 bvj[5]: 0.00000 svi[3]: 0.00000 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[0]: 0.211325 svi[4]: 0.211325 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[1]: 0.788675 svi[4]: 0.211325 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[2]: 0.00000 svi[4]: 0.211325 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[3]: 0.00000 svi[4]: 0.211325 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[4]: 0.00000 svi[4]: 0.211325 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[4]: 0.00000 bvj[5]: 0.00000 svi[4]: 0.211325 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[0]: 0.211325 svi[5]: 0.788675 svj[0]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[1]: 0.788675 svi[5]: 0.788675 svj[1]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[2]: 0.00000 svi[5]: 0.788675 svj[2]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[3]: 0.00000 svi[5]: 0.788675 svj[3]: 0.00000 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[4]: 0.00000 svi[5]: 0.788675 svj[4]: 0.211325 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::bvi[5]: 0.00000 bvj[5]: 0.00000 svi[5]: 0.788675 svj[5]: 0.788675 x[1]: 0.788675 0.00000 y[1]: 0.788675 0.00000 +DEAL::matrix = +DEAL::0.333333 0.166667 0.00000 0.00000 3.33333 1.66667 +DEAL::0.166667 0.333333 0.00000 0.00000 1.66667 3.33333 +DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL::0.00000 0.00000 0.00000 0.00000 0.00000 0.00000 +DEAL::33.3333 16.6667 0.00000 0.00000 333.333 166.667 +DEAL::16.6667 33.3333 0.00000 0.00000 166.667 333.333 -- 2.39.5