FECouplingValues.

[luca-heltai]

This is a draft PR that I would like to use as a discussion a class that can be used for general coupling operators.

I would like to compute bilinear forms of the following abstract type

$$ \int_{T_1} \int_{T_2} K(x_1, x_2) f(x_1) g(x_2) dT_1 dT_2, $$

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,
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:

$$ \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}, $$

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 two vectors 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:

$$ \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, $$

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 most common 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 provides 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 CouplingQuadratureType argument is usually chosen
to be CouplingQuadratureType::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 CouplingQuadratureType
argument (see the documentation of that enum class for more details).
An example usage of this class to assemble two coupled bilinear forms (for
example, for Boundary Element Methods) is the following:

... // double loop over cells that yields cell_1 and cell_2
fe_values_1.reinit(cell_1); fe_values_2.reinit(cell_2);
CouplingFEValues<dim> cfv(fe_values1, fe_values2,
                          CouplingQuadratureType::tensor_product);
FullMatrix<double> local_matrix (cfv.n_coupling_dofs,cfv.n_coupling_dofs);

// Extractor on left cell 
const auto v1 = cfv.left(FEValuesExtractor::Scalar(0)); 
const auto p1 = cfv.left(FEValuesExtractor::Scalar(1));
// Extractor on right cell 
const auto u2 = cfv.right(FEValuesExtractor::Scalar(0)); 
const auto q2 = cfv.right(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.coupling_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.coupling_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);
        }
    }
}

I would like some feedback on the interface before I dive into the implementation of more involuted complex things. The minimal test shows how the remapping of quadrature points and dof indices should be.

In principle, by augmenting the constructors of the current FEValuesViews with two renumbering vectors for dofs and quads (which are ignored if of zero lenght), I think we could use the full infrasctructure of the FEValuesViews which is already in place, without duplicating code.

Thoughts?

@pcafrica, @lwy5515, @fdrmrc

[luca-heltai]

@tjhei , would you mind giving a look at the snippet of code above and share your thoughts about the interface? Before I move forward and modify our FEvaluesViews classes, I'd like to have a consensus on how the interface should look like.

[tjhei]

I think this looks very reasonable. Is there a way to avoid the overhead of always having the permutation arrays present? Would it make sense to treat the identity differently. I expect this would cost performance otherwise.

[luca-heltai]

I'm not sure if this

const auto index = renumbering.empty() ? in_index :  renumbering[in_index];

is faster than unconditionally doing

const auto index = renumbering[in_index];

also for the identity. Opinions? I don't have a preference.

[luca-heltai]

So these are the steps needed for the concept presented here to work:

• Reorder an existing view (possibly without too many expensive operations) (see FEValuesViews::RenumberedView #15819)
• Create the correct dof reordering vectors
• Create the correct quadrature reordering vectors

[luca-heltai]

@tjhei , #15819 uses a different approach (inspired by your comment). There I do not touch the existing views, but I create a thin wrapper around them that encapsulates the renumbering. Let me know which approach you prefer.

[tjhei]

also for the identity. Opinions? I don't have a preference.

I would think a single if branch is the better choice over always storing the identity.

@tjhei , #15819 uses a different approach (inspired by your comment). There I do not touch the existing views, but I create a thin wrapper around them that encapsulates the renumbering. Let me know which approach you prefer.

Either option is fine with me. What do you prefer?

[luca-heltai]

Depends on #15819. Only the second commit introduces the actual FECouplingValues class.

[luca-heltai]

This is now fully functional.

• test 3: shows how to use this for BEM coupling (@lwy5515 )
• test 4: shows how to use this as a replacement for FEInterfaceValues (@fdrmrc)
• test 5: shows how to use this to compute bulk - surface coupling (@stefanozampini )

[fdrmrc]

The test fe_coupling_02.debug fails because the following assertion is not satisfied

dealii/include/deal.II/fe/fe_coupling_values.h

Lines 1515 to 1516 in 976d413

	Assert(dof_coupling_type != DoFCouplingType::independent,
	ExcMessage("Dofs are coupled. You cannot ask only for left dofs."));

How would you loop (within this interface) over DoF indices if DoFCouplingType::independent is selected?

[luca-heltai]

I added the assert later, and it was in the wrong direction. Now it should work.

[peterrum]

Looks good to me. Some minor comments. My biggest issue is that this might only work in serial (I am not counting p:s:T as parallel). The surface mesh and the bulk mesh are normally independent meshes, which are independently partitioned.

[peterrum]

These two classes look identical!? What is the reason?

[luca-heltai]

They are used to identify which of the two FEValuesBase you are addressing. Asking for fev[left] and fev[right] will have a different effect.

[luca-heltai]

Looks good to me. Some minor comments. My biggest issue is that this might only work in serial (I am not counting p:s:T as parallel). The surface mesh and the bulk mesh are normally independent meshes, which are independently partitioned.

This actually works perfectly fine on p::d::T (bulk) coupled with p::s::T (surface). Moreover, for BEM problems, once you set up the hierarchical matrices (based on shared trees), the p::s::T can be safely discarded. And one can create a surface p::f::T where the partitioning is induced by the (bulk) p::f::T (or p::d::T).

[kronbichler]

I think this looks good in principle, it is an exciting feature. However, I have two main concern, one regarding the naming left/right discussed further down, and one regarding the overall performance aspects. While I agree that the present class is nice, it is another incarnation (besides the FEValuesViews, FEInterfaceValues, etc) that fundamentally makes a slow-performance choice: For usability, we decide to return zeros whenever we are on the other side regarding the index. Do not get me wrong, this is not bad, but if one ever wants to us this in a context where performance matters, it seems this moves us in the wrong direction. This PR gets my approval because I do not have the resources to think for a good solution. Since we already have FEPointEvaluation or FEEvaluation in the matrix-free module that essentially provide 80% of what is required to assemble the coupling matrices, plus tensor product libraries like https://github.com/hyperdeal/hyperdeal that have algorithms for tensor product coupling, I strongly suggest that, in case performance starts to become relevant, to start looking there and talk to me.

[luca-heltai]

I think this looks good in principle, it is an exciting feature. However, I have two main concern, one regarding the naming left/right discussed further down, and one regarding the overall performance aspects. While I agree that the present class is nice, it is another incarnation (besides the FEValuesViews, FEInterfaceValues, etc) that fundamentally makes a slow-performance choice: For usability, we decide to return zeros whenever we are on the other side regarding the index. Do not get me wrong, this is not bad, but if one ever wants to us this in a context where performance matters, it seems this moves us in the wrong direction. This PR gets my approval because I do not have the resources to think for a good solution. Since we already have FEPointEvaluation or FEEvaluation in the matrix-free module that essentially provide 80% of what is required to assemble the coupling matrices, plus tensor product libraries like https://github.com/hyperdeal/hyperdeal that have algorithms for tensor product coupling, I strongly suggest that, in case performance starts to become relevant, to start looking there and talk to me.

I perfectly agree with your concerns. This is not thought to be fast for the computer, but fast for the programmer, similarly to FEInterfaceValues, and FEValuesViews. I wanted to have something with which it is easy to prototype bulk-surface, BEM, fractional Laplacian, etc.

Indeed, if we want to go fast, we need to take another route. However, I think that having this (slow) approach available is also necessary if one wants to be able to have a readable way to build these problems. I think this is a possible interface that follows the same steps we took elsewhere, and it would make it very easy to write coupled problems (with the same overhead of FEValuesViews/FEValues, etc.). If you see obvious ways in which this could be made faster, let's discuss it further. :)

[luca-heltai]

Ping @dealii/dealii ?

[bangerth]

Would you mind rebasing this to see whether the tests still all succeed?

[luca-heltai]

Only 68ab103 is relevant.

[bangerth]

Thanks. I will look at that over the next couple of days.

[bangerth]

Could you rebase now that the other patch has been merged? That should make the first commit go away.

[bangerth]

Part 1. The rest later today.

[bangerth]

Here's the rest. Nothing that will require much work.

I appreciate how extensively all of this is documented, and that you have 5 new tests. Good work.

[bangerth]

What would you think of putting these classes into a namespace FEValuesExtractors::FECouplingValues? This way they show up in the documentation in a similar place as the other extractors.

[luca-heltai]

@bangerth, I should have addressed all your comments. Thanks for your review!

[bangerth]

Thank you!