Galerkin BEM in deal.II

[Jakob Nielsen]

Dear deal.II community

I have attempted to implement a Galerkin BEM method by adapting the collocation method in Step-34. This first implementation is quite crude. I compute the singular integrals over identical elements i.e.
        int_T1 psi(x) int_T1 phi(y) G(x,y) dy dx
by wrapping the existing QGaussOneOverR in an outer loop over a QGauss quadrature. In all other cases I use a tensor QGauss quadrature. Although slow, the method does converges. However, I do not observe higher-order convergence when increasing the order of the basis functions. I assume this is due to the poor handling of the singular and near singular integrals.

I now want to improve this method focusing on the 3D problem with quadrilateral elements. The plan is to use one quadrature for the double integral itself through a series of coordinate transforms as proposed in [2, Sec 5.2]. Note that this is also the approach taken in BEM++ [1]. This method provides four different integration schemes (x_i,y_i,w_i), corresponding to the different cases

• T1 and T2 are identical.
• T1 and T2 share an edge, i.e. T1(s,0) = T2(s,0).
• T1 and T2 share a vertex, i.e. T1(0,0) = T2(0,0).
• T1 and T2 do not overlap.

In each case, the quadrature can be evaluated as
        sum_i w_i tilde(psi)(x_i) tilde(phi)(y_i) tilde(G)(x_i,y_i).

While reading the documentation, I came across FECouplingValues, which seems intended to support this type of element pair integration. However, I am struggling to understand how to practically handle the different geometric cases (identity, shared edge, shared vertex) and the required rotations within this framework.

Concretely, I would appreciate pointers on the following:

1. Is FECouplingValues the right tool for implementing these four quadrature cases for Galerkin BEM?
2. If so, is the intended workflow to explicitly construct paired (x_i,sqrt(w_i)) and (y_i,sqrt(w_i)) for each case (and possibly for each rotation), and then instantiate FECouplingValues separately for each configuration? This seems rather heavy, and I am wondering if there is a cleaner approach.
3. Do you know of any examples or code fragments in deal.II (or user projects) that implement Galerkin BEM?

References:
[1] BEM++ Documentation, http://www.bempp.org/quadrature.html
[2] Sauter and Schwab, Boundary Element Methods, 2011

Best regards,
Jakob

[Luca Heltai]

Dear Jakob,

your testcase is precisely one of the reasons why FECouplingValues was created. You should build a function that, given two cells, it returns the three vectors:

x, y, w

which you can then use as inputs to the FECouplingValues. Notice that FECouplingValues won’t help you with the reorientation of the cells/weights. You will have to do this yourself. FECouplingValues will help you later, once you have the points and weights.

Wenyu Lei had this implemented for triangles, if I recall correctly, using exactly the reference you pointed to, but I’m not finding the repository where this is committed now. Maybe Wenyu can chime in?

L.

> --
> The deal.II project is located at http://www.dealii.org/
> For mailing list/forum options, see https://groups.google.com/d/forum/dealii?hl=en
> ---
> You received this message because you are subscribed to the Google Groups "deal.II User Group" group.
> To unsubscribe from this group and stop receiving emails from it, send an email to dealii+un...@googlegroups.com.
> To view this discussion visit https://groups.google.com/d/msgid/dealii/a96009fc-7f70-4cba-adc6-dfbfe7168593n%40googlegroups.com.

[Jakob Nielsen]

Dear Luca,

Thank you for the quick and thoughtful response.

Based on your response, I have come up with the following skeleton for an assembly routine:

template <int dim> void GalerkinBEMProblem<dim>::assemble_system() {

for (const auto &cell_i : dof_handler.active_cell_iterators()) {

for (const auto &cell_j : dof_handler.active_cell_iterators()) {

const auto &[fev_i, fev_j] = fev_with_sauter_schwab_quadrature(cell_i, cell_j);

FECouplingValues<dim - 1, dim - 1, dim> cfev(

fev_i, fev_j, DoFCouplingType::independent,

QuadratureCouplingType::unrolled);

FullMatrix<Number> local_matrix = contribution_from_pair(cfev);

distribute_local_dofs_to_global(local_matrix, cell_i, cell_j);

}

}

};

In this case, fev_with_sauter_schwab_quadrature will handle all the logic with respect to choice and rotation of quadrature rule.

Is this how you suggest I utilise FECouplingValues?

One more thing:
Although the quadrature for the double integral can be written as (x_i,y_i,w_i) i=1,...,N, the x_i and y_i will be different combinations of a smaller subset of points. If we only included each point once, we could reduce the number of points by a factor of 1/8 (for identical cells) on each FEValues. However, this requires that the values are carefully rearranged when it comes to the stage of evaluating the quadrature. Perhaps support for this could be added through a new QuadratureCouplingType? Although I am not quite sure what the design should be, as we would also need a way to supply a method that specifies how the quadratures should be combined.

Best regards,
Jakob

[Luca Heltai]

Dear Jakob,

yes, this is roughly how we do it.

You have a list of possible quadrature point selection (and therefore fevalues), that can be combined together, and that you use to initialize the fevalues for cell_i and cell_j (that may well be the same).

I’m assuming this selection is done in a class that stores these guys (you don’t want to recreate them from scratch), and then returns with a member function (e.g., fev_with_sauter_schwab_quadrature) a reference to your two initialized FEValues.

However, in your case, you’d have to specify “contiguous”: both dof indices map to the same dofhandler, so you need to make sure that the loop over all indices is actually a loop over both cell_i and cell_j indices.

Your local matrix will look like this

A11 | A12
A =
A21 | A22

where A11 will contain interaction between dofs in cell_i, A22 those in cell j, and the off-diagonal will contain the interaction between dofs con cell_i and dofs on cell_j. The major role of FECouplingValues is to make this more readable, i.e.,

// Extractor on first cell
const auto extv = cfv.get_first_extractor(FEValuesExtractor::Scalar(0));

// Extractor on second cell
const auto extu = cfv.get_second_extractor(FEValuesExtractor::Scalar(0));

...
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[extv].value(i, q);
for (const unsigned int j : cfv.second_dof_indices())
{
const auto &u_j = cfv[extu].value(j, q);

local_matrix(i, j) += (Kernel(x_q, y_q) * v_i * u_j * cfv[0].JxW(q) * cfv[1].JxW(q));
}
}
}

fe_values_1.get_dof_indices(d1);
fe_values_2.get_dof_indices(d2);

auto coupling_dof_indices =
fecv.get_coupling_dof_indices(d1, d2);

constraints.distribute_local_to_global(local_matrix, coupling_dof_indices, global_matrix);



As far as your concerns about maybe saving a few quadrature points, the cost of this would be an additional indirection (i.e., you would have to map what unique quadrature point maps to what quadrature index). I’m unsure that saving in that direction would really be effective (remember: premature optimization is the source of all evils). Unless you know for a fact that it would be a problem to replicate a few quadrature points, I would not bother, at least in the first round of implementation. Once you have a working version, then you should profile it, and see where you are spending the most time.

It may be that this is going to be in the repeated quadrature points, but it may also be, in fact, that having those repeated but saving one indirection is faster.

Hope this helps.

L.

> To view this discussion visit https://groups.google.com/d/msgid/dealii/eaad257b-d3a2-4745-a556-834203d18911n%40googlegroups.com.