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.
> Il giorno 4 dic 2025, alle ore 10:26, Jakob Nielsen <[email protected]> ha
> scritto:
>
> 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
>
>
> tirsdag den 2. december 2025 kl. 16.34.18 UTC+1 skrev [email protected]:
> 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.
>
>
>
>
> > Il giorno 2 dic 2025, alle ore 15:04, Jakob Nielsen <[email protected]> ha
> > scritto:
> >
> > 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:
> > • Is FECouplingValues the right tool for implementing these four quadrature
> > cases for Galerkin BEM?
> > • 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.
> > • 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
> >
> > --
> > 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 [email protected].
> > To view this discussion visit
> > https://groups.google.com/d/msgid/dealii/a96009fc-7f70-4cba-adc6-dfbfe7168593n%40googlegroups.com.
> >
>
>
> --
> 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 [email protected].
> To view this discussion visit
> https://groups.google.com/d/msgid/dealii/eaad257b-d3a2-4745-a556-834203d18911n%40googlegroups.com.
--
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 [email protected].
To view this discussion visit
https://groups.google.com/d/msgid/dealii/33950411-424C-417E-80C2-A7FB6F4F929F%40gmail.com.
