Bembel: The Fast Isogeometric Boundary Element C++ Library for Laplace, Helmholtz, and Electric Wave Equation

In this article, we present Bembel, the C++ library featuring higher order isogeometric Galerkin boundary element methods for Laplace, Helmholtz, and Maxwell problems. Bembel is compatible with geometries from the Octave NURBS package and provides an interface to the Eigen template library for linear algebra operations. For computational efficiency, it applies an embedded fast multipole method tailored to the isogeometric analysis framework and a parallel matrix assembly based on OpenMP.


Introduction
The boundary element method (BEM) or method of moments (MoM) is a widely used tool for the solution of partial differential equations (PDEs) in engineering applications such as acoustic and electromagnetic scattering problems in homogeneous media. It is accepted in industry and implementations are available as software packages. Prominent examples are BEM++ [1] which is open-source and BETL [2] which is freely available for academic use. A commonality of these codes is that they are relying on mesh generators to create triangulation-based surface meshes consisting of flat triangles or lower order parametric elements. Thus, the order of convergence of (higher-order) boundary element methods is limited by the approximation error of the underlying mesh.
With the emergence of isogeometric analysis, see [3], boundary element methods have received increased attention by the community, since computer aided design (CAD) tools directly provide parametric representations of surfaces in terms of non-uniform rational B-splines (NURBS). This has recently led to a flourishing development of software concerning NURBS, see, e.g., [4,5]. Since the extraction of volume mappings from surface descriptions is an active research area with open problems, the use of isogeometric finite element methods is challenging in practice. Boundary element methods are the methods of choice in this setting, see also [6][7][8][9][10][11][12]. However, their implementation usually requires a significant amount of expert knowledge, which can lead non-experts to refrain from their usage.
The software library Bembel, Boundary Element Method Based Engineering Library, is a header-only library written in C++ (compliant to the 2011 standard) [13]. It provides a general framework and building blocks to solve boundary value problems governed by boundary integral operators within the isogeometric framework, for example the Laplace, Helmholtz, or electric wave equation. The development of the software started in the context of wavelet Galerkin methods on parametric surfaces, see [14], where the quadrature routines for the Green's function of the Laplacian have been developed and implemented. It was then extended to hierarchical matrices (H-matrices) in [15] and to H 2 -matrices and higher-order B-splines in [16]. With support of B-splines and NURBS for the geometry mappings, the Laplace and Helmholtz code became isogeometric in [7]. In [17], it has been extended to be applicable to non-scalar problems, by covering the electric field integral equation, and was modernized to a flexible C++ header-only library in early 2020.
The publication of this software package aims at making isogeometric boundary element methods available for a broader audience. Therefore, we aim at an easy to use C++ API, which streamlines the access to the underlying routines, with compatibility to the Eigen template library for linear algebra, see [18]. This is achieved, while still providing black-box H 2 -compression of the boundary element system matrices and OpenMP parallelization. The H 2 -compression yields an almost linear computational cost in the number of unknowns for assembly, storage requirements and matrix-vector multiplication of the system matrix. For the representation of geometries, Bembel features arbitrary parametric mappings, most prominently given as NURBS-mappings, which can directly be imported from files generated by the Octave NURBS package [4].
The structure of this document is as follows: Section 2 gives a brief introduction to isogeometric boundary elements, referring to the Appendix for details. Then, Section 3 explains our design considerations, after which Section 4 showcases a specific code example. With some familiarity to the important classes of our API, Section 5 and Section 6 then give an explanation of the structure of our library and its major building blocks. This is followed by a showcase of some simple numerical results in Section 7, obtained in complete analogy to the first code example. Afterwards, some remarks regarding the impact of our implementation are given in Section 8, and, finally, we conclude in Section 9.

Isogeometric boundary element methods
Bembel is able to treat general boundary value problems based on boundary integral operators, and provides examples to solve the well-known Laplace, Helmholtz, and electric wave equations as stated in Appendix A out of the box. For the Helmholtz and the electric wave equation, the provided examples assume nonresonant wave-numbers.
The boundary value problems can be solved by single layer potential ansatzes, see Appendix B, which require the solution of boundary integral equations by means of a numerical discretization. This and similar approaches are commonly known as boundary element methods. Bembel implements these through the use of a conforming Galerkin scheme.
One of the inherent advantages of boundary element methods over classical finite element techniques in the volume is that they can act directly on surface descriptions by NURBS from CAD programs, see Appendix C. This, in connection with corresponding spline spaces, leads to so-called isogeometric boundary element methods. Bembel implements these and assumes that the surface descriptions fulfil the requirements stated in Appendix D. The spline spaces for the Galerkin method are constructed as isogeometric multi-patch B-spline spaces, see Appendix E.
It is well known that in all provided examples solvability and uniqueness of the solution, both of the continuous problem as well as of its discrete counterpart, can be guaranteed by imposing reasonable assumptions on the boundary values, essentially guaranteeing ''physicality'' of the input data, see e.g. [19,20].
The computation of the system matrix requires the evaluation of singular integrals, see [20]. For the numerical quadrature of these integrals, we employ regularization techniques as described in [21]. The compression of the resulting densely populated system matrices is based on the embedded fast multipole method (FMM), which is tailored to the framework of isogeometric analysis, see [16,17], and fits into the framework of H 2 -matrices. Its particular advantage is that the matrix compression is directly applied on the reference geometry, that is, the unit square. Hence, the employed compression scheme profits from the inherently two-dimensional structure of the problem. The cost for assembly, storage requirement and matrix-vector multiplication for the system matrix are almost linear in the number of unknowns, see [16,17]. Moreover, this compression technique provably maintains the convergence behaviour for increasingly finer discretizations, cf. [16]. An in-depth mathematical analysis of the implemented approach together with numerical studies based on previous versions of Bembel is available in [7,16,17].

Design considerations
Most modern three-dimensional boundary element codes with built-in matrix compression are written in C or C++, resulting in efficient implementations. One of the central aims of Bembel is to provide a computationally efficient isogeometric boundary element code with a plain and simple user interface mimicking the mathematical setting. Therefore, the API of Bembel is designed in modern, template-based C++ and provides an interface for the Eigen template library for numerical linear algebra. This allows the user for a programming experience similar to Matlab and Octave and for the use of all matrix-free algorithms from the Eigen library. Particularly, all solvers provided by Eigen can be employed without further modifications.

Example program
Following the example of the Eigen library, Bembel is designed as a header-only library, structured into different modules. Before we discuss these modules in detail, we will present an example program to familiarize the reader with the typical program flow.
The discussed example is a main file that solves a Laplace problem and evaluates the potential at user defined points in space. The code corresponds to one of the examples provided in Bembel's repository, shortened by omitting loops over the polynomial degree and the refinement level, error measurements, and console output.
All of the presented functionality in the following example is either in the Eigen or Bembel namespace, thus lines 2 and 3 are merely for convenience. Line 4 initializes a Geometry object from a given file. As mentioned before, we support geometry files as exported by the NURBS package of Octave, and provide examples in the official repository that showcase how users can generate their own. MatrixXd gridpoints = // user defined points 7 8 std :: function < double ( Vector3d )> fun = 9 []( Vector3d in) { // user defined right hand side 10 }; 11 12 int polynomial_degree = // runtime parameter 13 int refinement_level = // runtime parameter 14 15 AnsatzSpace < LaplaceSingleLayerOperator > 16 ansatz_space (geometry , refinement_level , polynomial_degree ); 17 18 DiscreteLinearForm < DirichletTrace <double >, 19 LaplaceSingleLayerOperator > 20 disc_lf ( ansatz_space ); In line 15-16, the discrete space is initialized. The AnsatzSpace class generates the correct ansatz space for the LaplaceSingleLayerOperator, or possibly another linear operator defined by the user. The constructor of the AnsatzSpace class, offers a default parameter for knot repetition, which can be used for the refinement as well, thus making it possible to differentiate between p and k refinement.
Lines 18-22 assemble the right hand side of the linear system. The LinearForm, in this case given by the DirichletTrace implements a mapping that evaluates fun on the parametric surfaces w.r.t. the reference domain, which has to be passed to the LinearForm in line 21. Line 22 then calls the appropriate quadrature routines and assembles the right hand side.
Lines 24-26 take care of the assembly of the system matrix. In this example, we choose to assemble the system as our own and Eigen compatible H2Matrix format, but Eigen::MatrixXd may be used as well for a non-compressed assembly. The AnsatzSpace object is passed to the DiscreteOperator in line 25.
Therein, the compute() method then assembles the matrix. The specialization of the LinearOperatorBase class given via the LaplaceSingleLayerOperator needs to provide an evaluate_integrant_impl implementation, which is used to compute the matrix entries. A possible implementation is discussed below.
In lines 28-30, the linear system is solved by one of the matrix-free solvers provided by the Eigen library.
Afterwards, in lines 32-35, a DiscretePotential object is created, which evaluates the pointwise solution to the PDE from the computed Cauchy data. Therein, all required information is passed to the classes via the specialization LaplaceSingleLayerPotential, which is once again user defined.

Modules and notes on their use
Now, that a basic familiarity with a typical program flow has been established, we will introduce all Bembel modules in alphabetical order and discuss their purpose.

AnsatzSpace
The AnsatzSpace module contains the routines managing the discrete space on the surface of the geometry. Specifically, this is realized through the four classes SuperSpace, Projector, Glue, and AnsatzSpace. Therein, SuperSpace manages local polynomial bases on every element. Through a transformation matrix generated by the template class Projector which depends on the specialization of LinearOperatorBase and its defined traits, the SuperSpace can be related to B-Spline bases on every patch. To build conforming spaces (in the case of DifferentialForm::DivergenceConforming through continuity of the normal component across patch interfaces, in the case of DifferentialForm::Continuous through global C 0 -continuity), the template class Glue assembles another transformation matrix to identify degrees of freedom across edges. Then, a coefficient vector in the SuperSpace can be related to one of the smooth B-Spline bases, as explained in [9,17].

ClusterTree
The ClusterTree module introduces a cluster hierarchy on parametric surfaces. The contained ClusterTree class takes a Geometry and introduces an element structure on it, in the form of an ElementTree. Each ElementTreeNode of the ElementTree corresponds to an entire parametric mapping (at the trees first level) or to a sub element induced by a recursive refinement strategy. The leaves of the tree correspond to the elements on which the SuperSpace introduces shape functions.

DuffyTrick
The DuffyTrick module provides quadrature routines for (nearly) singular integrals. Therein, the implementation is as described in the appendix of [14], essentially being an adaptation of the so-called Sauter-Schwab quadrature rules [21]. These routines, together with the Quadrature module, can easily be used as a standalone module on any type of quadrilateral discretization.

DummyOperator
This module provides a specialization of LinearOperatorBase for testing.

Geometry
The Geometry module provides the functionality required for a parametric geometry description in a computational framework. The Geometry class is the interface between geometry description and the remainder of the code. We provide an implementation of NURBS discretized patches via the Patch class, but the code is easily extensible to other parametric descriptions mapping from [0, 1] 2 to parts of the geometry, as long as the corresponding methods for point evaluation and evaluation of the pointwise Jacobian are implemented. The Geometry module depends only on the Spline module and can be used standalone for the implementation of custom numerical codes.

H2Matrix
The H2Matrix module provides functionality for an efficient compression of the system matrix and reduction of the computational cost of the matrix-vector multiplication. Therein, the implemented algorithms correspond to the ones presented in [7,16].

Helmholtz
The Helmholtz module provides specializations to solve Helmholtz problems.

IO
The IO module provides input-output functionality, including routines for VTK file export, timing, and writing log files.

Laplace
The Laplace module provides specializations to solve Laplace problems.

LinearForm
The LinearForm template class must be specialized for the assembly of the right hand side. Exemplarily, Bembel provides an implementation of the scalar DirichletTrace and the vectorvalued RotatedTangentialTrace required to solve Laplace and Helmholtz Dirichlet problems and the electric field integral equation. Other specializations implementing different approaches can be easily implemented by appropriate modifications within the corresponding header files.

LinearOperator
Specializations of the LinearOperatorBase class govern the behaviour of most other modules. To provide a valid specialization, methods for evaluation of the integrand, i.e., including the test functions, must be provided. Moreover, the corresponding specialization of LinearOperatorTraits must be provided, allowing other classes to determine crucial properties such as the numerical type of the problem (in general double or std::complex<double >) and the type of discretization, i.e., either DifferentialForm:: Continuous, corresponding to a discrete subspace of H 1/2 , DifferentialForm::DivConforming, corresponding to a discrete subspace of H H H −1/2 × (div Γ , Γ ), or DifferentialForm::Discontinuous, corresponding to a discrete subspace of H −1/2 (Γ ).

Maxwell
The Maxwell module example specializations required to solve problems based on the electric wave equation.

Potential
The Potential module introduces means to evaluate the solution to a PDE via a suitable integral operator taking the unknown of the linear system as input. It relies, once again, on a suitable specialization corresponding of the LinearOperatorBase class.

Quadrature
The Quadrature module provides quadrature routines for [0, 1] n for n = 1, 2. The implementation is for arbitrary quadrature orders, where for higher-orders the required data is generated by the solution of the corresponding three term recurrence. The Quadrature module can be used standalone.

Spline
The Spline module provides basic routines related to spline function and local polynomials. The polynomials are implemented through template recursion, and should be interfaced through the SuperSpace class if used within Bembel. The Spline module is independent of other modules and can easily be used as the basis of other numerical codes.

User specific implementation
In all of the above, the user needs to provide certain implementations in order to specify custom operators. These are the routines for the evaluation of the integrand, as well as the definition of certain traits. We will briefly discuss the implementation of the Laplace single layer case, following the previous code example. 1 class LaplaceSingleLayerOperator ; This is an example of LinearOperatorTraits. They pose a way for core algorithms to deduce types, such as complex or real values, and matrix formats associated with the problem to be solves. The OperatorOrder enum enables the core routines to choose an appropriate quadrature degree, and the Form designates that the basis is to be assembled as a discretization of H −1/2 (Γ ), i.e., discontinuous across patch boundaries. Further traits designate information for the H2Matrix compression, for a discussion of these, we refer to the Doxygen documentation of the code.
For an implementation of the matrix assembly, so-called SurfacePoints are passed as an input. A SurfacePoint p can be handled as follows. The SurfacePoint is a typedef for an Eigen::Matrix<double,Eigen:: Dynamic,12>, where the first two components s hold the coordinate in the reference element, the third component w an optional (quadrature) weight, the fourth to sixth entry f the corresponding point on the geometry in R 3 , and the remaining components f_dx and f_dy the tangential vectors corresponding to differentiation in x and y direction. A surface point can be ''unwrapped'' via a code snipped as above, generating references to the segments of the vector. Now, an implementation of evaluateIntegrand_impl could look as follows. (* intval ) += super_space . basisInteraction ( s, t) * evaluateKernel (x_f , y_f)* x_kappa * y_kappa * ws * wt; 21 22 return ; 23 } Lines 1-11 show that this is a specialization of the LinearOperatorBase class. In the implementation of the evaluation of the integrand, lines 13-19 merely extract information from the given input, the SurfacePoints and the SuperSpace. The SuperSpace is aware of the refinement structure and yields methods for basis function evaluation. The actual implementation relevant for matrix assembly is in line 20, which directly corresponds to an evaluation of the matrix local element matrix for quadrature points and weights encoded in the SurfacePoints.
Note that the implementation of a corresponding LaplaceSingleLayerPotential is completely analogous, and a specialization for the FMM is provided as an example in the repository as well, within the Laplace module.

Numerical results
Running a more elaborate version of the example code discussed above, including for loops of the polynomial order and the refinement level, as well as an error measurement via manufactured solution, yields the numerical results illustrated in Fig. 1, which behave as predicted by theory, cf. [7]. We remark that the visualization of the density in Fig. 1 was generated by the routines provided in the IO module. More involved numerical examples computed with previous versions of Bembel, including Helmholtz and Maxwell problems, are discussed in [7,17].

Impact
The implementation of boundary element methods in three spatial dimensions is a non-trivial task due to the necessary numerical evaluation of singular integrals and the required matrix compression to achieve computational efficiency. While this is already true for lowest-order implementations for the Laplace equation on boundary triangulations with flat elements, difficulties increase on isogeometric (or parametric) surfaces, isogeometric B-spline (or higher-order) boundary element spaces, and more involved electromagnetic problems. These implementationrelated issues lead many people to refrain from the use of boundary element methods, even when a specific engineering problem is known to be solved best therewith. Our software package aims to provide a state-of-the-art boundary element toolbox for engineers who would like to apply competitive isogeometric boundary element methods in a black-box fashion. This allows users to freely employ boundary element methods as a tool in involved engineering applications. To the best of our knowledge, the functionality of higher-order B-spline boundary element spaces is the first open source implementation available for three dimensions.
The major strengths of Bembel are the capability of implementing general boundary integral operators, and the direct integration into the Eigen Linear Algebra Library. This allows the user for a straightforward pre or post processing of data with a clean user interface. The dense and matrix-free algorithms of Eigen provide different kinds of solvers for the Galerkin systems or eigenvalue problems.

Conclusion
Bembel is an open-source library enabling users to apply isogeometric boundary element methods for potential, acoustic, electromagnetic, and many other problems in a black-box fashion. This is achieved through an easy-to-use API both for black-box use, as well as implementation of custom operators, and compatibility with the Eigen template library for linear algebra. Moreover, it incorporates state-of-the-art compression techniques for large problems, as well as OpenMP parallelization.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.

Acknowledgments
This work is supported by DFG Grants SCHO1562/3-1 and KU1553/4-1 within the project Simulation of superconducting cavities with isogeometric boundary elements (IGA-BEM and the Helmholtz equation with Sommerfeld radiation conditions towards infinity. In both cases, the boundary values have to be understood in the usual sense of traces, see, e.g. [20]. In addition, Bembel can treat the electric wave equation with Silver-Müller radiation conditions towards infinity, which can be derived from Maxwell's equations. Again, the boundary values have to be understood in the sense of traces, see [19].

Appendix B. Boundary integral equations
The three boundary value problems (LP), (HP), and (MP) can each be solved by means of a single layer potential ansatz, i.e., setting and (MS). We note that S L and S H are isomorphisms from H −1/2 (Γ ) to H 1/2 (Γ ), whereas S M is an isomorphism on H −1/2 × (div Γ , Γ ), see [19,20] for details. Thus, in view of a conforming Galerkin method, piecewise polynomial boundary element spaces are sufficient for the boundary element based solution of (LP) and (HP), whereas (MP) requires divergence conforming boundary element spaces. Unique solvability of the corresponding linear systems can be proven, see [19,20]. towards infinity, see [19,20].