𝖧𝖫𝖨𝖡𝗉𝗋𝗈 implements various tools for discretization and solving boundary element methods, namely grid generation, function spaces and bilinear forms.

The implementation is based on solving the integral equation

\begin{equation*} \int_{\Gamma} k(x,y) \mathbf{u}(y) dy = \mathbf{f}(x), \quad x \in \Gamma \end{equation*}

with \(\Gamma = \partial \Omega \subset \mathbf{R}^3\) and \(\mathbf{u} : \Gamma \to \mathbf{R}\) being sought for a given right hand side \(\mathbf{f} : \Gamma \to \mathbf{R}\).

The solution should be computed using Galerkin discretisation with ansatz functions \(V = \{\phi_i, 0 \le i < n\} \) and test functions \(W = \{\psi_i, 0 \le i < m\} \), which results in a linear equation system

\[A u = f\]

where \(u\) contains the coefficients of the discretised \(\mathbf{u}\) and \(A\) is defined by

\[ a_{ij} = \int_\Gamma \int_\Gamma \phi_i(x) k(x,y) \psi_j(y) dy dx \]

The right hand side \(f\) is given by

\[ f_i = \int_{\Gamma} \psi_i(x) f(x) dx. \]

Remarks: Please note, that the functionality of the implementation is limited, e.g. only triangular grids, piecewise constant and linear function spaces. Normally this sub library is used for benchmarking H-matrix algorithms.

Please see Boundary Element Methods for a full example (Helmholtz).

Grids

Grids can be loaded from files. Here different file formats are supported, e.g., the HLIB, PLY, GMSH and SurfaceMesh format, with the restriction of triangular grids. One can either use file format auto detection and use the general function

std::unique_ptr< TGrid >

read_grid ( const std::string & filename );

or directly the read function from the classes TAutoGridIO, THLibGridIO, TPlyGridIO, TSurfMeshGridIO and TGMSHGridIO.

Furthermore, basic grids can be constructed with one of

auto make_sphere () -> std::unique_ptr< TRefinableGrid >;
auto make_cube   () -> std::unique_ptr< TRefinableGrid >;
auto make_square () -> std::unique_ptr< TRefinableGrid >;

which will construct a sphere (diameter 2, centered at \((0,0,0)^T\)), cube ( \((0,0,0)^T \ldots (1,1,1)^T\)) or square ( \((0,0,0)^T \ldots (1,1,0)^T\)) grid.

As the return type TRefinableGrid suggest, these grids can be refined using the function

std::unique_ptr< TRefinableGrid > refine () const;

This can be used to refine the grid until the number of triangles is sufficiently large:

auto  grid = make_sphere();
 
for ( uint  i = 0; i < nlevels; ++i )
    grid = grid->refine();

Remarks: Refinement is not supported with grids loaded from files!

Function Spaces

Based on these grids functions spaces are constructed, defining the degrees of freedom (DoF). 𝖧𝖫𝖨𝖡𝗉𝗋𝗈 implements TConstFnSpace for piecewise constant functions, TLinearFnSpace for piecewise linear functions and TConstEdgeFnSpace for constant normal linear tangential (CN/LT) edge elements (for Maxwell applications).

During construction, the single argument is the grid:

TConstFnSpace cfnspace( grid.get() );

TLinearFnSpace lfnspace( grid.get() );

As function spaces define the DoFs, their number is available with the n_indices function

std::cout << cfnspace.n_indices() << std::endl;

std::cout << lfnspace.n_indices() << std::endl;

Furthermore, the coordinates of all DoFs are provided by the function spaces using the build_coord function:

auto ccoord = cfnspace->build_coord();

auto lcoord = lfnspace->build_coord();

which can be used for clustering.

Bilinear Forms

𝖧𝖫𝖨𝖡𝗉𝗋𝗈 implements the bilinear forms for the following BEM kernels:

Laplace single \(\left(\frac{1}{\|x-y\|_2}\right)\) and double \(\left( \frac{\langle n, y-x \rangle}{\|x-y\|_2^3} \right)\) layer potential (TLaplaceSLPBF and TLaplaceDLPBF)
Helmholz single \(\left( \frac{e^{i\kappa \|x-y\|_2}}{\|x-y\|_2}\right)\) and double \(\left( \frac{e^{i \cdot \kappa \|x-y\|_2} (i\cdot \kappa \|x-y\|_2 - 1) \langle n, y-x \rangle}{\|x-y\|_2^3} \right)\) layer potential (THelmholtzSLPBF and THelmholtzDLPBF)
Maxwell electric and magnetic field integral equation (TMaxwellEFIEBF and TMaxwellMFIEBF)
Mass matrix (TMassBF, TMaxwellEFIEMassBF and TMaxwellMFIEMassBF)
exponential kernel \(\left(e^{ - \|x-y\|_2}\right)\) (TExpBF)

The implementation of all these is based on quadrature (Sauter-Schwab triangular quadrature), the order of which is am optional parameter for the corresponding constructors. Beside this, the functions spaces for the ansatz and test function spaces are the other parameters:

TLaplaceSLPBF< TConstFnSpace, TLinearFnSpace > laplace_bf( & cfnspace, & lfnspace );

Right Hand Side

The computation

\[ f_i = \int_{\Gamma} \psi_i(x) f(x) dx. \]

of the right hand side is also implemented using quadrature and implemented for the scalar valued types in TQuadBEMRHS.

For this, the user needs to implement the function \(f\) by deriving from TBEMFunction, which defines an eval function for a point \(x\) and normal direction \(n\):

virtual value_t eval ( const T3Point & x,

const T3Point & n ) const;

This function is abstract in TBEMFunction and needs to be implemented in derived classes:

class TMyRHS : public TBEMFunction< real_t >
{
public:
    virtual real_t
    eval  ( const T3Point &  x,
            const T3Point &  n ) const
    {
        return ... ;
    }
};

With the help of TQuadBEMRHS, the discrete version of \(f\) can be computed:

TMyRHS                                rhsfn;
TQuadBEMRHS< TConstFnSpace, real_t >  rhs_build( quad_order );
auto                                  rhs = rhs_build.build( & cfnspace, & rhsfn );

RHS for Maxwell

For Maxwell computations, the RHS is vector valued. A special type complex3_t is available for holding three complex values. With this, the user defined RHS function looks like

class TMyRHS : public TBEMFunction< complex3_t >
{
public:
    virtual complex3_t
    eval  ( const T3Point &  x,
            const T3Point &  n ) const
    {
        return complex3_t{ ..., 
                           ...,
                           ... };
    }
};

The computation of the RHS vector is performed by TMaxwellEFIERHS and TMaxwellMFIERHS.