At its core, the assembly of boundary integral operators in Bempp is based on an efficient evaluation of integrals of the form

I = \int_{T_1}\int_{T_2}g(\mathbf{x},\mathbf{y})\phi(\mathbf{y})\overline{\psi(\mathbf{x})}\mathrm{d}\mathbf{y}\mathrm{d}\mathbf{x},

where \displaystyle g(\mathbf{x},\mathbf{y}) is a weakly singular kernel, \displaystyle T_1 and \displaystyle T_2 are triangles, and \displaystyle \psi and \displaystyle \phi are basis functions on \displaystyle T_1 and \displaystyle T_2. We need to distinguish two cases:

  • \displaystyle T_1 and \displaystyle T_2 are fully disjoint,
  • \displaystyle T_1 and \displaystyle T_2 share a vertex, edge or are the same triangle.

In the first case a fully regular Gauss quadrature rule on triangles can be used, that is we approximate the integral as

I\approx \sum_{i=1}^{N_1}\sum_{j=1}^{N_2}g(\mathbf{x}_i,\mathbf{y}_j)\phi(\mathbf{y}_j)\overline{\psi(\mathbf{x}_i)}\omega^{(1)}_i\omega^{(2)}_j,

where \mathbf{x}_i\in T_1 and \mathbf{y}_j\in T_2 are quadrature points with weights \omega^{(1)}_i and \omega^{(2)}_j.

In the second case the singularity needs to be taken into account. This is done using Duffy type transformations. Details of the singular integration can be found in the book Boundary Element Methods by Sauter and Schwab.

How many integration points are used?

The cost of the assembly of boundary element matrices is dominated by the number of kernel evaluations for each pair of triangles \displaystyle (T_1,T_2) in the approximation of the corresponding Galerkin integral.

For disjoint triangles a standard tensor Gauss rule based on symmetric Gauss points over triangles is used. The following table shows the order, the number of points used on a triangle and the total number of kernel evaluations over the pair \displaystyle (T_1, T_2).

Order Gauss points in triangle Total number of kernel evaluations
1 1 1
2 3 9
3 4 16
4 6 36
5 7 49
6 12 144
7 13 169
8 16 256
9 19 361
10 25 625
11 27 729
12 33 1089
13 37 1369
14 42 1764
15 48 2304
16 52 2704
17 61 3721
18 70 4900
19 73 5329
20 79 6241

If the triangles \displaystyle T_1 and \displaystyle T_2 are not disjoint the total number of kernel evaluations depends on how the two triangles intersect each other. Let \displaystyle O be the singular integration order. Then the total number of kernel evaluations \displaystyle E_K is given as follows.

  • The triangles are adjacent at a single vertex: \displaystyle E_K = 2 \left\lfloor (O+2)/2\right\rfloor^4
  • The triangles are adjacent at an edge: \displaystyle E_K = 5 \left\lfloor (O+2)/2\right\rfloor^4
  • The triangles coincide (\displaystyle T_1==T_2) \displaystyle E_K = 6 \left\lfloor (O+2)/2\right\rfloor^4

The default singular quadrature order in Bempp is 6. Hence, the total number of kernel evaluations for an element pair is 512 in the case of adjacent vertices, 1280 in the case of adjacent edges and 1536 in the case of coincident triangles. This is significantly more than in the case of disjoint triangles. However, the number of non-disjoint triangle pairs depends only linearly on the mesh size and this typically only takes a small fraction of the overall computing time.

Controlling the quadrature order in Bempp

Bempp provides fine-grained control over the quadrature order. The simplest way to change the quadrature order is to use the global parameter object bempp.api.global_parameters. For example, to change the singular quadrature order to 5, use the following command.

bempp.api.global_parameters.quadrature.double_singular = 5

The quadrature order for approximating the regular integral over pairs of disjoint triangles \displaystyle T_1 and \displaystyle T_2 is dependent on the normalised distance between them. The normalised distance is defined as the distance between the centers of the triangles divided by the maximum of the areas of \displaystyle T_1 and \displaystyle T_2.

Three integration zones are defined, near, medium and far. The following parameters are set as default.

Zone Max normalised distance Double Order Single Order
near 2 4 4
medium 4 3 3
far infinity 2 2

Each quadrature rule is used for pairs of triangles with normalised distance up to the rule’s max normalised distance. The single order is the quadrature order used for the evaluation of potentials, where the normalised distance calculated between the evaluation point and the center of the triangle. For the discretisation of grid functions, the far single order is used. The values given here are relatively sane defaults. However, we note that good values for the quadrature orders depend on the geometry and the accuracy requirements and may be very different from the default values stated here.

To change for example the values of the medium distance quadrature orders the following commands can be used.

bempp.api.global_parameters.quadrature.medium.max_rel_dist = 5
bempp.api.global_parameters.quadrature.medium.single_order = 1
bempp.api.global_parameters.quadrature.medium.double_order = 1