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

TL;DR

Bembel is a C++ library that implements high-order isogeometric boundary element methods for solving Laplace, Helmholtz, and Maxwell problems efficiently using fast multipole and parallel assembly techniques.

Contribution

It introduces a fast, efficient C++ library integrating isogeometric boundary element methods with advanced computational techniques.

Findings

Achieves high computational efficiency for boundary element problems

Supports complex geometries via compatibility with NURBS

Utilizes fast multipole method for scalable performance

Abstract

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.

Tables

Table 1. Table 1: Code metadata

Nr.	Code metadata description	Please fill in this column
C1	Current code version	v0.9
C2	Permanent link to code/repository used for this code version	http://www.bembel.eu, cf. [13]
C3	Legal Code License	GPL3
C4	Code versioning system used	Git
C5	Software code languages, tools, and services used	C, C++, OpenMP
C6	Compilation requirements, operating environments & dependencies	CMake, C++11, Eigen Linear Algebra Library
C7	If available Link to developer documentation/manual	http://www.bembel.eu/ $\to$ Documentation
C8	Support email for questions	info@bembel.eu

Equations

$$\begin{array}{rlrl}-\Delta u_{\mathrm{L}}=& 0& & \text{in}\Omega ,\\ u_{\mathrm{L}}=& g_{\mathrm{L}}& & \text{on}\Gamma ,\end{array}$$

$$\begin{array}{rlrl}-\Delta u_{\mathrm{H}}-{\kappa }_{\mathrm{H}}^2{u}_{\mathrm{H}}=& 0& & \text{in }{\Omega }^c,\\ u_{\mathrm{H}}=& g_{\mathrm{H}}& & \text{on }\Gamma ,\end{array}$$

$$\begin{array}{rlrl}\mathrm{curl}\mathrm{curl}{\mathbf{E}}_{\mathrm{M}}-{\kappa }_{\mathrm{M}}^2{\mathbf{E}}_{\mathrm{M}}=& 0& & \text{in}{\Omega }^c,\\ \mathbf{n}\times {\mathbf{E}}_{\mathrm{M}}=& {\mathbf{g}}_{\mathrm{M}}& & \text{on }\Gamma ,\end{array}$$

$$u_{\mathrm{L}}(\mathbf{x})=$$

$$u_{\mathrm{H}}(\mathbf{x})=$$

$${\mathbf{E}}_{\mathrm{M}}(\mathbf{x})=$$

$$+\frac{1}{{\kappa }_{\mathrm{M}}^2}{\nabla }_{\mathbf{x}}{\int }_{\Gamma }\frac{e^{-i{\kappa }_{\mathrm{M}}\Vert \mathbf{x}-\mathbf{y}{\Vert }_2}}{4\pi \Vert \mathbf{x}-\mathbf{y}{\Vert }_2}{\mathrm{div}}_{\Gamma }{\mathbf{j}}_{\mathrm{M}}(\mathbf{y})d{\sigma }_{\mathbf{y}},\mathbf{x}\in {\Omega }^c.$$

$${\mathcal{S}}_{\mathrm{L}}{\rho }_{\mathrm{L}}=g_{\mathrm{L}},{\mathcal{S}}_{\mathrm{H}}{\rho }_{\mathrm{H}}=g_{\mathrm{H}},{\mathcal{S}}_{\mathrm{M}}{\rho }_{\mathrm{M}}={\mathbf{g}}_{\mathrm{M}}$$

$$b_j^p(x)$$

$${\mathbf{\gamma }}_j(x,y):=\sum _{0\le j_1<k_1}\sum _{0\le j_2<k_2}\frac{{\mathbf{c}}_{j_1,j_2}{b}_{j_1}^{p_1}(x)b_{j_2}^{p_2}(y)w_{j_1,j_2}}{{\sum }_{i_1=0}^{k_1-1}{\sum }_{i_2=0}^{k_2-1}{b}_{i_1}^{p_1}(x)b_{i_2}^{p_2}(y)w_{i_1,i_2}},$$

$${\mathbf{\gamma }}_i:\square \to {\Gamma }_i\text{ with }{\Gamma }_i={\mathbf{\gamma }}_i(\square )\text{ for }i=1,2,\dots ,M,$$

$$\Gamma =\bigcup _{i=1}^M{\Gamma }_i,$$

$${\mathbb{S}}_{\mathbf{p},\mathbf{\Xi }}^0(\square )$$

$${\mathbb{S}}_{\mathbf{p},\mathbf{\Xi }}^1(\square )$$

$${\mathbb{S}}_{\mathbf{p},\mathbf{\Xi }}^2(\square )$$

$${\mathbb{S}}_{\mathbf{p},\mathbf{\Xi }}^0(\Gamma )\subset H^{1/2}(\Gamma ),{\mathbb{S}}_{\mathbf{p},\mathbf{\Xi }}^1(\Gamma )\subset {\mathbf{H}}_{\times }^{-1/2}({\mathrm{div}}_{\Gamma },\Gamma ),{\mathbb{S}}_{\mathbf{p},\mathbf{\Xi }}^2(\Gamma )\subset H^{-1/2}(\Gamma ).$$

Full text

††journal: arXiv

1 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 accessible 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 unsolved open problems, the use of isogeometric finite element methods is quite 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 written in C and C++ [13]. It solves boundary value problems governed by the Laplace, Helmholtz or electric wave equation within the isogeometric framework. The development of the software started in the context of wavelet Galerkin methods on parametric surfaces, see [14], where the integration routines for the Green’s function of the Laplacian have been developed and implemented. It was then extended to hierarchical matrices ($\mathcal{H}$-matrices) in [15] and to $\mathcal{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]. Finally, in [17], it has been extended to the electric field integral equation.

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 $\mathcal{H}^{2}$-compression of the boundary element system matrices and OpenMP parallelized matrix assembly. The $\mathcal{H}^{2}$-compression yields an almost linear complexities in the number of unknowns for assembly, storage requirements and matrix-vector multiplication for the system matrix. For the representation of geometries, Bembel features arbitrary parametric mappings, most prominently given as NURBS-mappings and can directly be imported from files generated by the Octave NURBS package [4].

The structure of this document is as follows. Section 2 shortly recalls the main concepts of isogeometric boundary element methods. Section 3 is concerned with design considerations for the library, whereas Section 4 deals with the implementation. Afterwards, in Section 5, we discuss an example program, whereas Section 6 discusses the significance of the code. Finally, in Section 7, we will conclude our review.

2 Isogeometric Boundary Element Methods

Bembel is able to treat potential, acoustic scattering, and electromagnetic scattering problems governed by the well-known Laplace, Helmholtz, and electric wave equations as stated in A. For the Helmholtz and the electric wave equation, Bembel assumes non-resonant wave-numbers.

The boundary value problems can be solved by single layer potential ansatzes, see 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 C. This, in connection with corresponding spline spaces, leads to so-called isogeometric boundary element methods. Bembel implements this and assumes that these surface descriptions fulfill the requirements stated in D. The spline spaces for the Galerkin method are constructed as isogeometric multi-patch B-spline spaces, see E.

It is well known that in all cases solvability and uniquenes 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 $\mathcal{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 complexity 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].

3 Design Considerations

Most modern three-dimensional boundary element codes with built-in matrix compression are written in C or C++, which are proven to lead to 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. Therefore, the API of Bembel is designed in modern 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 iterative solvers provided by Eigen can be employed without further modifications. In contrast, the low-level routines of Bembel are written in optimized C code which is tuned towards efficiency. Shared memory parallelism by the use of OpenMP is provided and hence aims at multicore machines.

4 Implementation

4.1 Geometry

With Bembel, geometry representations can directly be imported from files generated by the Octave NURBS package [4]. However, we should mention that the implementation of the Geometry class is rather flexible such that the NURBS mappings may be replaced by arbitrary smooth parametric mappings defined by the user.

4.2 PDE

The PDE to solve and the potential ansatz to use are defined by a derivation of the base class PDEproblem. The currently available derivations are LaplaceSingle, HelmholtzSingle, and MaxwellSingle, which implement the boundary integral equations (BIE) for Laplace, Helmholtz and Maxwell problems listed in B. The derived classes contain all necessary information for the PDE and the boundary integral equation under consideration, see A and B. Note that, for HelmholtzSingle and MaxwellSingle, the construction of a class instance requires a wave number described as std::complex<double>.

4.3 Discretization

Given a Geometry object and a class inheriting from PDEproblem, e.g. LaplaceSingle, an object of the template class Discretization<LaplaceSingle> describes a conforming boundary element space on refinement level L and polynomial degree P. The type of the B-spline boundary element space is automatically defined through the class inheriting from PDEproblem, e.g. LaplaceSingle.

4.4 Boundary Values

The load vectors of the Galerkin systems can be obtained via quadrature from a std::function object describing the boundary values. For the Laplace case, computeRhs takes a Discretization and a std::function<double(Eigen::Vector3d)>) object and returns the load vector as an Eigen::VectorXd. For the Helmholtz and Maxwell case, the return value is an Eigen::VectorXcd. Note that in this case the std::function object has to be modified accordingly, i.e. it must accept arguments of type Eigen::Vector3cd. Bembel provides exemplary functions. They can be wrapped into this framework, e.g, via C++ lambdas.

4.5 System Matrix

Given a specialized Discretization, e.g.  Discretization<LaplaceSingle>, and a polynomial degree for the $\mathcal{H}^{2}$-compression, a corresponding specialization of the template class HierarchicalMatrix, e.g.  HierarchicalMatrix<LaplaceSingle>, may be used to assemble the compressed system matrix of the Galerkin system. Note that the type of the PDE and the geometry are directly encoded into the discretization, such that no further information is required here.

The template class HierarchicalMatrix<T> inherits from the corresponding specialization of the EigenBase class, i.e. EigenBase<HierarchicalMatrix<T>>. Therefore, a HierarchicalMatrix object may use the Eigen interface and can be passed to all functions which accept EigenBase<Derived> objects as parameters.

The Eigen methods rows() and cols() are specialized in order to return the dimensions. Moreover, to facilitate the multiplication of Eigen::VectorX objects, the specialization Product<HierarchicalMatrix, Rhs, AliasFreeProduct> of the Eigen Product class is provided. In particular, the operator * is overloaded to support expressions like H * x, where H is a HierarchicalMatrix<PDEproblem> object, and x an Eigen::VectorX object, with a fitting Scalar type.

4.6 Solution of the System of Linear Equations

The solution of linear systems within Bembel can be performed by the matrix-free solvers from Eigen. To that end, the template class HierarchicalMatrixinherits the Eigen::internal::traits of an Eigen::SparseMatrix with the same Scalar type as in the corresponding example from the Eigen documentation on matrix-free solvers. We opted for the Eigen::SparseMatrix traits, since they seem to impose the least requirements. Hence, in view of the implemented general matrix-vector (GEMV) product H * x, all Eigen iterative solvers which solely rely on this operation may be used.

4.7 Evaluation of the Solution

Having the density vector rho from the iterative solver, the solution to the PDE in either the exterior or the interior domain can be evaluated by the function Sol::evalSolution(gridpoints, rho, myDisc). In here, gridpoints should be an Eigen::Matrix<double, Eigen::Dynamic, 3> matrix of n three-dimensional points in the domain of interest. Depending on the PDE, sol is either an Eigen::VectorXd, Eigen::VectorXcd or Eigen::Matrix<std::complex<double>, Eigen::Dynamic, 3> object, where each row contains the solution at the corresponding gridpoint.

4.8 Visualization

The Bembel::Vis::plotDiscretizationToVTK functions can be used to visualize rho and the geometry in vtk-format. Their functionality is explained extensively in the src/examples/example_Visualization.cpp file.

5 Example Program

We discuss a full example for the solution of an electromagnetic scattering problem, which is located in src/examples/example_MaxwellSingle.cpp.

⬇

1int main() {

2 using namespace Bembel;

3 using namespace Eigen;

4

5 const int multipoleDegree = 16;

6 const int knotRepetition = 1;

7 const std::complex<double> wavenumber(3, 0);

8

9 Geometry geom("geometry.dat");

10 MaxwellSingle pde(wavenumber);

11 MatrixXd gridpoints = Util::makeSphereGrid(6, 10);

12

13 const std::function<Vector3cd(Vector3d, std::complex<double>)>

14 fun = ... ; // Define incident wave

15

16 for (int P = 1 ; P < 10 ; ++P)) {

17 for (int L = 0 ; L < 6 ; ++L) {

18 Discretization<MaxwellSingle> disc(geom, pde, P, knotRepetition, L);

19 VectorXcd rhs = Rhs::computeRhs(disc, fun);

20 HierarchicalMatrix<MaxwellSingle> hmat(disc, multipoleDegree);

21 GMRES<HierarchicalMatrix<MaxwellSingle>, IdentityPreconditioner> gmres;

22 gmres.compute(hmat);

23 VectorXcd rho = gmres.solve(rhs);

24 MatrixXcd pot = Sol::evalSolution(gridpoints, rho, disc);

25 }

26 }

27 return (0);

28}

At the beginning of the program, the polynomial degree for the $\mathcal{H}^{2}$-compression, the knot repetition of the B-spline spaces, and the wave number are defined. The Geometry object is loaded from the file "geometry.dat", and the PDE is defined for a given wave number. Util::makeSphereGrid(6,10) generates ten evaluation points on ten lines of latitude for the scattered field on a sphere around the origin with radius 6. Afterwards, an incident wave for the scattering problem is defined in fun.

The rest of the example program loops over different polynomial degrees P and refinement levels L. In the loop, the discrete space is defined as a Discretization object, the load vector and the system matrix are assembled and solved with the GMRES solver from the Eigen (unsupported) package. Finally, the scattered field is evaluated at the previously defined points.

Note that the program is straightforwardly adapted to the Laplace and Helmholtz cases as indicated by the example files in the src/examples folder. These examples also provide error measurements of the solution to the underlying PDE and output of these errors into log-files. These were not included to the above example for a better understanding of the main functionality of the code. LaTeX-files creating error plots such as illustrated in Figure 1 can be found in the LaTeX folder. The routines can visualize the intermediate quantity rho as illustrated in Figure 1.

6 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 implementation-related 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.

A major strongpoint of Bembel is 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 matrix-free algorithms of Eigen provide different kinds of iterative solvers for the Galerkin systems or eigenvalue problems.

7 Conclusion

Bembel is an open-source library enabling users to apply isogometric 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 and compatibility with the Eigen template library for linear algebra. Moreover, it encorporates state-of-the-art compression techniques for large problems, as well as OpenMP parallelization.

Acknowledgement

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). The work of Jürgen Dölz was partially supported by SNSF Grants 156101 and 174987, as well as the Graduate School of Computational Engineering at TU Darmstadt and the Excellence Initiative of the German Federal and State Governments and the Graduate School of Computational Engineering at TU Darmstadt. Michael Multerer was partially supported by SNSF Grant 137669. The work of Felix Wolf is supported by the Excellence Initiative of the German Federal and State Governments and the Graduate School of Computational Engineering at TU Darmstadt.