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 a quadrature rule

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}}\langle w,{V}_{L}u\rangle ={\displaystyle \underset{\Gamma}{\int}}{\displaystyle \underset{\Gamma}{\int}}w\left(x\right){U}_{L}(y-x)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}[k,\ell ]=\langle {w}_{k},{V}_{L}{u}_{\ell}\rangle ={\displaystyle \underset{\mathrm{supp}\left({w}_{k}\right)}{\int}}{\displaystyle \underset{\mathrm{supp}\left({u}_{\ell}\right)}{\int}}{K}_{k\ell}({U}_{L},{w}_{k},{u}_{\ell})\phantom{\rule{0.278em}{0ex}}d{s}_{y}d{s}_{x}\phantom{\rule{0.278em}{0ex}},\phantom{\rule{2.00em}{0ex}}{K}_{k\ell}({U}_{L},{w}_{k},{u}_{\ell})={w}_{k}\left(x\right){U}_{L}(y-x){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/galerkinquadrature.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}{\left|z\right|}\phantom{\rule{0.167em}{0ex}},\phantom{\rule{2.00em}{0ex}}{U}_{H}(z;\kappa )=\genfrac{}{}{0.1ex}{}{1}{4\pi}\genfrac{}{}{0.1ex}{}{\mathrm{exp}(-\kappa |z\left|\right)}{\left|z\right|}\phantom{\rule{0.167em}{0ex}},\phantom{\rule{2.00em}{0ex}}z\in {\mathbb{R}}^{3}\phantom{\rule{0.278em}{0ex}}.$

Note that the wave-number
$\kappa \in \u2102$
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: * type of tau_y, regular, coincident, edge-adjacent, vertex-adjacent rules * type of tau_x, regular, coincident, edge-adjacent, vertex-adjacent rules */ 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

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}[k,\ell ]=\langle {\mathbf{w}}_{k},{V}_{H}{\mathbf{u}}_{\ell}\rangle ={\displaystyle \underset{\mathrm{supp}\left({\mathbf{w}}_{k}\right)}{\int}}{\displaystyle \underset{\mathrm{supp}\left({\mathbf{u}}_{\ell}\right)}{\int}}{K}_{k\ell}({U}_{H},{\mathbf{w}}_{k},{\mathbf{u}}_{\ell})\phantom{\rule{0.278em}{0ex}}d{s}_{y}d{s}_{x}\phantom{\rule{0.278em}{0ex}},\phantom{\rule{2.00em}{0ex}}{K}_{k\ell}({U}_{H},{\mathbf{w}}_{k},{\mathbf{u}}_{\ell})={\mathbf{w}}_{k}\left(x\right)\cdot \left({U}_{H}\right(y-x\left){\mathbf{u}}_{\ell}\right(y\left)\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`

`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 typedef GalerkinQuadrature< element_t, element_t, 7, 7, 25, 25, 64, 64, // sledgehammer 9, 9 > quadrature_t;

A calculation with the recompiled code reveals the following result for the Laplace equation's single layer potential

1.8546e-01 7.4551e-02 7.4551e-02 7.4551e-02 7.4551e-02 7.9824e-02 3.9251e-02 3.9251e-02 7.4551e-02 3.9251e-02 7.9824e-02 3.9251e-02 7.4551e-02 3.9251e-02 3.9251e-02 7.9824e-02

Now, the matrix looks fine. Nevertheless, the number of 64 Gaussian points is ridiculous. However, this example already illustrates the sensitivity of the integration's accuracy with respect to the chosen quadrature rule. We will comment on the integration routines in more detail in one of the upcoming examples.

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.