At its core, the assembly of boundary integral operators in Bempp is based on an efficient evaluation of integrals of the form
where is a weakly singular kernel, and are triangles, and and are basis functions on and . We need to distinguish two cases:
- and are fully disjoint,
- and 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
where and are quadrature points with weights and .
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 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 .
|Order||Gauss points in triangle||Total number of kernel evaluations|
If the triangles and are not disjoint the total number of kernel evaluations depends on how the two triangles intersect each other. Let be the singular integration order. Then the total number of kernel evaluations is given as follows.
- The triangles are adjacent at a single vertex:
- The triangles are adjacent at an edge:
- The triangles coincide ()
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.
The quadrature order for approximating the regular integral over pairs of disjoint triangles and 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 and .
Three integration zones are defined, near, medium and far. The following parameters are set as default.
|Zone||Max normalised distance||Double Order||Single Order|
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.single_order = 1
bempp.api.global_parameters.quadrature.medium.double_order = 1