## Example 3: Generating Boundary Layer potentials

Welcome to Boundary Element Methods! Contrary to the former examples which dealt more or less with the setup of the discretisation model this example is the very first one which is devoted to the discretisation of Boundary Integral Operators. Again, we will utilise the former example from the section called “Example 2: The DofHandler concept” to build an application which computes two discrete boundary layer potentials as they have been theoretically introduced in the section called “Trace Operators and Boundary Integral Operators”. Thereby, the discretised bilinear forms will correspond to those finite element spaces which have been introduced previously. The complete source code for this example may be found at `\$BETL_ROOT/tutorial/example_3/main.cpp`.

In detail, the current example deals with

• Defining fundamental solutions

• Defining kernel functions

• Defining integrators

• Defining BEM operators

For example the bilinear corresponding to the Laplace equation's single layer potential reads as

Equation 4.1.

$\phantom{\rule{1.00em}{0ex}}〈w,{V}_{L}u〉=\underset{\Gamma }{\int }\underset{\Gamma }{\int }w\left(x\right){U}_{L}\left(y-x\right)u\left(y\right)\phantom{\rule{0.278em}{0ex}}d{s}_{y}d{s}_{x}$

where we have skipped the trace operators for simplicity. The Galerkin discretisation of this bilinear form results in the matrix

Equation 4.2.

$\phantom{\rule{1.00em}{0ex}}{\mathsf{V}}_{L}\left[k,\ell \right]=〈{w}_{k},{V}_{L}{u}_{\ell }〉=\underset{\mathrm{supp}\left({w}_{k}\right)}{\int }\underset{\mathrm{supp}\left({u}_{\ell }\right)}{\int }{K}_{k\ell }\left({U}_{L},{w}_{k},{u}_{\ell }\right)\phantom{\rule{0.278em}{0ex}}d{s}_{y}d{s}_{x}\phantom{\rule{0.278em}{0ex}},\phantom{\rule{2.00em}{0ex}}{K}_{k\ell }\left({U}_{L},{w}_{k},{u}_{\ell }\right)={w}_{k}\left(x\right){U}_{L}\left(y-x\right){u}_{\ell }\left(y\right)\phantom{\rule{0.278em}{0ex}}.$

The computation of matrices in form of ${\mathsf{V}}_{L}$ is what BETL is mainly about. To achieve this in a user convenient fashion, it tries to mimic the mathematical notation by respective data structures. For instance, the expression ${K}_{k\ell }$ in the equation above is represented by the data structure ``` GalerkinKernel```. As we will see, the declaration of this data type follows directly its mathematical definition.

To begin with, here are the additional headers which have to be included:

```#include "traits/complextraits.hpp"
#include "fundsol/fundsol.hpp"
#include "kernel/galerkinkernel.hpp"
#include "integration/galerkinintegrator.hpp"
#include "bem_operator/bem_operator.hpp"
```

Following Equation 4.2 we can conclude that we need to define the fundamental solution ${U}_{L}$. We do not need to define one of the functions ${w}_{k}$ or ${u}_{\ell }$ since these functions already have been determined by the `FEBasis` and the `DoFHandler`. In BETL the fundamental solution is defined via two enumerators. While the first enumerator determines the underlying partial differential operator the second one defines which boundary potential layer should be instantiated. The following code

```// declare two fundamental solution data types
typedef FundSol< LAPLACE  , SLP > fs_lapl_slp_t;
typedef FundSol< HELMHOLTZ, SLP > fs_helm_slp_t;
```

declares two fundamental solutions. The first one is nothing but the already-known fundamental solution for the Laplace equation's single layer potential, the second one corresponds to the single layer potential for the Helmholtz operator. Their definitions are given below

Equation 4.3.

$\phantom{\rule{1.00em}{0ex}}{U}_{L}\left(z\right)=\genfrac{}{}{0.1ex}{}{1}{4\pi }\genfrac{}{}{0.1ex}{}{1}{|z|}\phantom{\rule{0.167em}{0ex}},\phantom{\rule{2.00em}{0ex}}{U}_{H}\left(z;\kappa \right)=\genfrac{}{}{0.1ex}{}{1}{4\pi }\genfrac{}{}{0.1ex}{}{\mathrm{exp}\left(-\kappa |z|\right)}{|z|}\phantom{\rule{0.167em}{0ex}},\phantom{\rule{2.00em}{0ex}}z\in {ℝ}^{3}\phantom{\rule{0.278em}{0ex}}.$

Note that the wave-number $\kappa \in ℂ$ is a complex number by what the evaluation result of the Helmholtz fundamental solution will be a complex number. Later in this example we will comment on how BETL treats complex numbers.

Now we can go ahead with the definition of kernel. In BETL the kernel is nothing but the concatenation of the test- and trial-functions together with the fundamental solution. Here is the declaration of two kernels

```// declare laplace kernel based on the Lagrangian based fe space
typedef GalerkinKernel< fs_lapl_slp_t,
feb_e3_const_disc_lagr_t,
feb_e3_const_disc_lagr_t > kernel_lapl_t;
// declare helmholtz kernel based on the edge based fe space
typedef GalerkinKernel< fs_helm_slp_t,
feb_e3_lin_const_div_t,
feb_e3_lin_const_div_t > kernel_helm_t;
```

Clearly, the declarations of the kernel data types follow directly the kernel definition from Equation 4.2.

In general, an analytic integration of the kernel function in Equation 4.2 is impossible. Hence, in most instances we are forced to replace that integration by an appropriate quadrature rule. Contrary to the classical finite element schemes, the quadrature in boundary element methods is by no means straightforward. Due to the singularities of the kernel functions we need to distinct between these four different cases

• the distance of all points on two elements is in any case positive (regular case)

• two elements are identical (coincident case)

• two elements share a common edge (edge adjacent case)

• two elements share a common point (vertex adjacent case).

A general quadrature which tackles those four different cases can be, e.g, found in [StrSchwb11]. The quadrature scheme presented there covers weak and strong singularities for triangular as well as for quadrilateral elements – whether these elements are curved or not. Hence, it represents the Swiss army knife for the evaluation of kernel functions in 3-dimensional Galerkin-based boundary element methods. Within BETL this special quadrature is the the default integration scheme and it is declared by the following code

```// declare a quadrature rule
/*
* template paramters:
*/
typedef GalerkinQuadrature< element_t, 7, 25, 16, 9
element_t, 7, 25, 16, 9 > quadrature_t;
```

The `GalerkinQuadrature` data type takes up to 10 template arguments. The first four template arguments specify the quadrature rule for the inner integration on ${\tau }_{y}$. Beside the element type one has to prescribe the quadrature rules for the regular, coincident, the edge-adjacent, and for the vertex-adjacent quadrature. The remaining four template parameters specify the quadrature for the outer integration on ${\tau }_{x}$. Typically, these quadrature rules equal those for the innter integration. Therefore, the last four template parameters are optional and can be omitted. For the regular integration on triangular elements BETL uses quadrature formulas as they are, e.g, given in [Dnvnt85]. Refer to `gausstria.hpp` for available quadrature rules on triangles. Independent of the element type, the quadrature for the singular integrations is based on tensor Gauss point rules. Therefore, the remaining pairs of integer numbers in `GalerkinQuadrature` need to be square numbers defining the Gaussian points on a quadrilateral reference element. The implemented tensor product rules are specified in the header `gausstensor.hpp`.

With a specific quadrature rule and a particular kernel type we are ready to declare the integration data type named `GalerkinIntegrator`. These declarations are quite simple and read as

```// declare the integrator types
typedef GalerkinIntegrator< kernel_lapl_t, quadrature_t > integrator_lapl_t;
typedef GalerkinIntegrator< kernel_helm_t, quadrature_t > integrator_helm_t;
```

Now, it remains to define instances of some of the above declared data types. The definition of the Helmholtz fundamental solution in Equation 4.3 features the wave-number as a parameter. As already mentioned this parameter is in general a complex number. Since C++ offers no native complex data type there might be more than just one complex data type implementation around. Therefore, BETL offers access to complex data types via `ComplexTraits`. Note that the complex data type is chosen by the `ComplexTraits` with help of precompiler directives.

```// define a wave number
typedef ComplexTraits< double >::complex_type complex_t;
const complex_t kappa( 0., 1. );
```

Now, what follows are the definitions of the fundamental solution, of the kernel, and of the integrator. Thereby, the former instance serves as a constructor argument for the latter one. The following source code illustrates this.

```// create fundsol, kernel, and integrator instances
fs_lapl_slp_t     fs_lapl_slp;
fs_helm_slp_t     fs_helm_slp( kappa );
kernel_lapl_t     kernel_lapl( fs_lapl_slp );
kernel_helm_t     kernel_helm( fs_helm_slp );
integrator_lapl_t integrator_lapl( kernel_lapl );
integrator_helm_t integrator_helm( kernel_helm );
```

With the instantiation of the integrators we have finished the setup of the more BEM specific tasks. We are now able to perform the integration on two particular boundary elements ${\tau }_{x}$ and ${\tau }_{y}$. In other words, the integrator is a functor which is responsible for the computation of one or of a couple of specific matrix entries. Since we aim at the computation of the complete system matrix ${\mathsf{V}}_{L}$ we need a structure which combines the `DoFHandler` (which stores the information on how to assemble the degrees of freedom) with the `GalerkinIntegrator` (which actually is able to evaluate the matrix entries). Within BETL this structure is called `BemOperator` and its declaration as well as its instantiation is

```// declare the discrete bem operators
typedef BemOperator< integrator_lapl_t,
dofhandler_lagr_t,
dofhandler_lagr_t > bem_operator_lapl_t;
typedef BemOperator< integrator_helm_t,
dofhandler_div_t,
dofhandler_div_t > bem_operator_helm_t;
// instantiate the discrete bem operators
bem_operator_lapl_t bem_operator_lapl( integrator_lapl,
dofhandler_lagr, dofhandler_lagr );
bem_operator_helm_t bem_operator_helm( integrator_helm,
dofhandler_div, dofhandler_div );
```

Finally, the boundary layer potentials are generated via the `compute()` method

```// compute the boundary potentials
bem_operator_lapl.compute( );
bem_operator_helm.compute( );
```

While the command `bem_operator_lapl.compute( )` computes the single layer potential for the Laplace equation as it is stated in Equation 4.2 the command `bem_operator_helm.compute( )` evaluates the following bilinear form

Equation 4.4.

$\phantom{\rule{1.00em}{0ex}}{\mathsf{V}}_{H}\left[k,\ell \right]=〈{\mathbf{w}}_{k},{V}_{H}{\mathbf{u}}_{\ell }〉=\underset{\mathrm{supp}\left({\mathbf{w}}_{k}\right)}{\int }\underset{\mathrm{supp}\left({\mathbf{u}}_{\ell }\right)}{\int }{K}_{k\ell }\left({U}_{H},{\mathbf{w}}_{k},{\mathbf{u}}_{\ell }\right)\phantom{\rule{0.278em}{0ex}}d{s}_{y}d{s}_{x}\phantom{\rule{0.278em}{0ex}},\phantom{\rule{2.00em}{0ex}}{K}_{k\ell }\left({U}_{H},{\mathbf{w}}_{k},{\mathbf{u}}_{\ell }\right)={\mathbf{w}}_{k}\left(x\right)\cdot \left({U}_{H}\left(y-x\right){\mathbf{u}}_{\ell }\left(y\right)\right)\phantom{\rule{0.278em}{0ex}}.$

Note that the test- and trial-space in Equation 4.4 are based on Equation 2.3.

The generations of the discrete boundary layer potentials ${\mathsf{V}}_{L}$ and ${\mathsf{V}}_{H}$ almost conclude this section. But as in the previous section, for debug purposes it is necessary to allow the boundary potentials to be exported either directly to the screen or to a file. The way the output of those potentials is done equals the output of the dof handler objects in the previous section. We simply apply the stream operator to the underlying matrix structures of the respective bem operators. Here is the code for getting the references to the generated matrices.

```// get references to the created matrix structures
bem_operator_lapl_t::const_reference V_lapl = bem_operator_lapl.giveMatrix( );
bem_operator_helm_t::const_reference V_helm = bem_operator_helm.giveMatrix( );
```

With the references `V_lapl` and `V_helm` the output of the potentials can be done in the following way:

```// write matrices to files stream
const std::string fname_lapl = basename + "_LAPL.dat";
const std::string fname_helm = basename + "_HELM.dat";
std::ofstream out_lapl( fname_lapl.c_str() );
std::ofstream out_helm( fname_helm.c_str() );
out_lapl << V_lapl;
out_helm << V_helm;
out_lapl.close( );
out_helm.close( );
```

Finally, building this example via `make betl_ex3` and calling `betl_ex3 tetrahedron_4.msh` produces the two output files `tetrahedron_4_LAPL.dat` and `tetrahedron_4_HELM.dat`, respectively. Those files can then be, e.g., directly processed with GNU Octave. The additional on-screen output for this example shows the constructor messages of the `BemOperator` and some statistics for its `compute()` method.

```<<<-------------------------------------------------------------------
BemOperator
----------------------------------------------------------------------
dimension:                4 x 4
symmetry:                 symmetric
(dense) mem. consumption: 7.62939e-05 MB
precision:                Double precision
numerical type:           d
acceleration method:      NO_ACCELERATION
------------------------------------------------------------------->>>

<<<-------------------------------------------------------------------
BemOperator
----------------------------------------------------------------------
dimension:                6 x 6
symmetry:                 symmetric
(dense) mem. consumption: 0.000320435 MB
precision:                Double precision
numerical type:           4compIdE
acceleration method:      NO_ACCELERATION
------------------------------------------------------------------->>>

<<<-------------------------------------------------------------------
BemOperator::compute()
----------------------------------------------------------------------
computational time:     0.001672 sec
performed integrations: 16
saved integrations:     0
------------------------------------------------------------------->>>

<<<-------------------------------------------------------------------
BemOperator::compute()
----------------------------------------------------------------------
computational time:     0.006681 sec
performed integrations: 16
saved integrations:     0
------------------------------------------------------------------->>>
```

The output may differ on your system since the numerical type is queried via the `typeid` function from the `typeinfo` header. Moreover, the computational times will surely also differ on your system.

Now, we will have a closer look on the computational results. The entries of ${\mathsf{V}}_{L}$ are stored in the file `tetrahedron_4_LAPL.dat` with its content given below.

```   1.8546e-01   7.4552e-02   7.4552e-02   7.4552e-02
7.4542e-02   7.9824e-02   3.9244e-02   3.9250e-02
7.4542e-02   3.9250e-02   7.9824e-02   3.9244e-02
7.4542e-02   3.9244e-02   3.9250e-02   7.9824e-02
```

Since the bilinear form in Equation 4.1 is symmetric and since we apply a symmetric Galerkin discretisation to it we would assume the system matrix to be symmetric as well. Unfortunately, this is not the case in this example. As one already might assume, this is due to the singular integrations. The mesh contains four elements which build a tetrahedron. Thus, the integration routines which are called are only the singular integration schemes for the coincident and edge adjacent cases. In this example the coincident integrations are only called to compute the diagonal entries, all other entries are evaluated by the edge adjacent integration routines. And these edge adjacent integration schemes rely on unsymmetric quadrature schemes. In order to allay this phenomenon we will increase the number of Gaussian points. In doing so we crack a nut with a sledgehammer and prescribe the following quadrature.

```// declare a quadrature rule
7, 7,
25, 25,
64, 64,   // sledgehammer
9,  9
```   1.8546e-01   7.4551e-02   7.4551e-02   7.4551e-02
One final comment has to be given on the fact that obviously all entries of the matrix ${\mathsf{V}}_{L}$ have been evaluated and that no advantage has been taken out of the bilinear form's symmetry properties. Boundary element methods rely on non-local operators resulting in fully populated system matrices. Both, in terms of memory consumption as well as in terms of computational costs this leads to a quadratic complexity of the algorithm. For practical applications this complexity cannot be tolerated. Thus, we need to incorporate the so-called fast boundary element methods in order to tackle real world problems. From this point of view the dense matrices which have been created here can be considered only as some kind of debug matrices. For instance, they may serve to check some of the system's properties or they may be used to validate the fast boundary element methods which will be introduced later. Therefore, it simply does not make sense to take advantage of the bilinear form's symmetry at this stage.