Very brief BEM4I tutorial

Installation

The current version is only available as a header-only library. There is no need for compilation and installation, you only have to place the project's directory into a convenient location and include required headers into your code. However, BEM4I requires Eigen and METIS external libraries for its functionality. Therefore, after downloading BEM4I, download and install the libraries and provide the correct paths during the compilation (see the example makefiles for Intel and GCC compilers in the examples directory).

Usage

You can find example code inside the examples directory. Use the EXAMPLE variable in the examples/main.cpp file to select one of the problems from the provided list (currently for the Laplace and Helmohltz equations). Simple makefiles are provided for the Intel and GCC compilers in the same directory.

Let us briefly describe the workflow using a simple code for the Helmholtz Dirichlet problem located in the examples/Helmholtz/testHelmholtzDirichlet.h source file. To compile the example, use #define EXAMPLE 200 in the main.cpp file. Structure of the remaining examples is similar.

Solution of the Dirichlet problem for the Helmohltz equation

We intend to solve

where $a=-\sqrt{3}\kappa, b=2\kappa$, and $\kappa=2$.

Within the main function we call the testHelmholtzDirichlet function with the following arguments

• input surface mesh file,
• number of refinements by quadrisection,
• number of refinements by dividing each side to thirds,
• bool indicating whether to map the mesh to unit ball,
• type of quadrature - 0 for semi-analytical, 1 for fully numerical,
• quadrature order,
• number of points for evaluation by the representation formula.

One of the first things done in the testHelmholtzDirichlet is creating the instance of the SurfaceMesh3D class and loading the input mesh using the load method. BEM4I currently supports triangular surface meshes stored in text files with the following structure:

      3
      3

      8
      1.0 -1.0 -1.0
      1.0 1.0 -1.0
      ...

      12
      0 1 5
      0 5 4
      

The first two lines denotes the spatial dimension and number of nodes per element, respectively. The following number denotes the total number of nodes in the mesh and is followed by the coordinates of individual nodes. Finally, number of elements and list of elements (indices into the node list) follows. The indexing starts with 0. The elements nodes $\mathbf{x}_1, \mathbf{x}_2, \mathbf{x}_3$ are listed such that $(\mathbf{x}_2-\mathbf{x}_1)\times(\mathbf{x}_3-\mathbf{x}_1)$ defines the outer normal to the element.

After loading, the mesh is refined and possibly mapped to the unit sphere. Next, an array of quadrature orders with size two for semi-analytical approach (see, e.g., Rjasanow, Steinbach. The Fast Solution of Boundary Integral Equations, Springer 2007) or four for fully-numerical approach (see, e.g., Sauter, Schwab. Boundary Element Methods, Springer 2004) is allocated. Quadrature orders for integration over disjoint elements are set using the int quadDisj [] = { 4, 4 }; line. The values in the array denotes quadrature orders for outer and inner integrals, respectively.

Next, discrete counterparts of boundary element spaces are created over the loaded mesh:

      BESpace< LO, SC > bespace00( &mesh, p0, p0 );
      BESpace< LO, SC > bespace10( &mesh, p1, p0 );
      BESpace< LO, SC > bespace11( &mesh, p1, p1 );
      

After the mesh argument, the constructor takes the type of ansatz and test functions, respectively. Currently, piece-wise constant (p0) and piece-wise linear (p1) basis functions are supported. Here we therefore create space bespace00 used for construction of the single layer matrix, space bespace10 used for construction of the double layer matrix, and space bespace11 used for construction of the discretized identity operator.

Next, we create empty objects of the FullMatrix class and discrete bilinear forms that will take care of the assembly of the matrices:

      FullMatrix< LO, SC > * V = new FullMatrix< LO, SC >( 0, 0 );
      BEBilinearFormHelmholtz1Layer< LO, SC > formV( &bespace00, quadOrder, kappa, quadType, quadDisj );

      FullMatrix< LO, SC > * K = new FullMatrix< LO, SC >( 0, 0 );
      BEBilinearFormHelmholtz2Layer< LO, SC > formK( &bespace10, quadOrder, kappa, quadType, quadDisj );
      

The constructors of the bilinear forms take the discrete boundary spaces, quadrature orders for close elements, wavenumber, quadrature type and orders of disjoint quadrature, respectively, as arguments.

To assemble the required matrix, one only has to call the assemble method of the appropriate bilinear form object.

      formK.assemble(*K);
      

In the next section of the code we integrate over the boundary of the computational domain and compute the right hand side using

      IdentityOperator< LO, SC > id( &bespace10 );
      id.apply( dir, rhs, false, 0.5, 0.0 );
      K->apply( dir, rhs, false, 1.0, 1.0 );
      

The A.apply(x, y, trans, alpha, beta) method of the IdentityOperator and FullMatrix classes computes the matrix-vector multiplication in the form z = beta*y + alpha*A*x. After assembling the matrix $V$ we therefore obtain the system in the form

and solve it using the GMRES algorithm

      V->GMRESSolve( rhs, sol, 1e-12, 1000, 1000 );
      

Besides the right-hand-side and solution vectors, here we set the relative precision, maximum number of iterations, and number of iterations before restart.

In the remaining part of the code we check the relative error of the solution, evaluate the solution inside of the computational domain using the representation formula, and store the result on the boundary as well as the result in a point-cloud in the ParaView file format.

Licence

Copyright (c) 2021, VSB - Technical University of Ostrava All rights reserved.

Redistribution and use in source and binary forms, with or without modification, are permitted provided that the following conditions are met:

• Redistributions of source code must retain the above copyright notice, this list of conditions and the following disclaimer.
• Redistributions in binary form must reproduce the above copyright notice, this list of conditions and the following disclaimer in the documentation and/or other materials provided with the distribution.
• Neither the names of VSB - Technical University of Ostrava nor the names of its contributors may be used to endorse or promote products derived from this software without specific prior written permission.

THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS “AS IS” AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL VSB - TECHNICAL UNIVERSITY OF OSTRAVA BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.