HLIBpro  2.9.1
Boundary Element Methods

In this example, an integral equation is to be solved by representing the discretised operator with an H-matrix, whereby the H-matrix is constructed using a matrix coefficient function for entries of the equation system. Furthermore, for the iterative solver a preconditioner is computed using H-LU factorisation.

# Problem Description

Given 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 Galerkin discretisation with ansatz functions $$V = \{\phi_i, 0 \le i < n\}$$ and test functions $$W = \{\psi_i, 0 \le i < m\}$$ leads to 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.$

# Surface Grid and Function Spaces

The code starts with the standard initialisation:

#include <iostream>
#include <hlib.hh>
using namespace HLIB;
using namespace std;
int main ( int argc, char ** argv )
{
try
{
INIT();

๐ง๐ซ๐จ๐ก๐๐๐ implements grid generation for sphericial and cubical grids. Furthermore, it supports triangular surface grids stored in several file formats, e.g. HLIB, PLY, SurfaceMesh and Gmsh v2 format. Although individual I/O classes for each file format exist, you may also use automatic file format detection:

auto grid = make_grid( "grid.tri" );

Based upon the grid, function spaces for the ansatz and the test space are defined. For Laplace and Helmholtz kernels, piecewise constant and linear function spaces are available, whereas for Maxwell kernels, piecewise constant edge space (RWG elements) is implemented.

auto ansatz_fnspace = make_unique< TLinearFnSpace >( grid.get() );
auto test_fnspace = make_unique< TLinearFnSpace >( grid.get() );
Remarks
Different grids may be used for ansatz and test spaces.

The functions spaces provide necessary geometrical information for construction cluster trees and block cluster trees for the defined index sets:

auto coord1 = ansatz_fnspace->build_coord();
auto coord2 = test_fnspace->build_coord();
TAutoBSPPartStrat part_strat;
TBSPCTBuilder ct_builder( & part_strat );
auto ct1 = ct_builder.build( coord1.get() );
auto ct2 = ct_builder.build( coord2.get() );
complex kappa( -2, 1 );
TBCBuilder bct_builder;
auto bct = bct_builder.build( ct1.get(), ct2.get(), & adm_cond );

Here, the admissibility condition THiLoFreqGeomAdmCond for oscillatory kernels was used, which tests if the given number of wave lengths (10) will fit into a given cluster to enable low-rank approximations. This is especially necessary if kappa is large.

# Definition of Kernel Function and Matrix Construction

In ๐ง๐ซ๐จ๐ก๐๐๐, different kernels are defined by special bilinear forms, each derived from TBEMBF. For reasons of efficiency, e.g. for basis function evaluation, the ansatz spaces are provided to each bilinear form as template arguments. For the Helmholtz single layer potential, the corresponding bilinear form is declared as:

using slpbf_t = THelmholtzSLPBF< TLinearFnSpace, TLinearFnSpace >;
slpbf_t slpbf( kappa, ansatz_fnspace.get(), test_fnspace.get() );

Here, kappa is the wave number of the underlying problem.

Remarks
When template based classes are used, as here with the Helmholtz bilinear form, type aliases will increase readability of the code.

For matrix construction, the bilinear form is not fully sufficient as actual matrix coefficients are needed. These are provided by the coefficient function TBFCoeffFn using the bilinear form together with TPermCoeffFn to convert the different orderings:

TBFCoeffFn< slpbf_t > slp_coeff_fn( & slpbf );
TPermCoeffnFn< complex > coeff_fn( & slp_coeff_fn, bct->row_ct()->perm_i2e(), bct->col_ct()->perm_i2e() );

Finally, the low-rank approximation technique and the block-wise accuracy gave to be defined, being ACA+ and $$\epsilon = 10^{-4}$$ in our case. Equipped with these, the TDenseMBuilder class can construct the discretised Helmholtz single layer potential:

TACAPlus< complex > aca( & coeff_fn );
TTruncAcc acc( 1e-4, 0.0 );
TDenseMBuilder< complex > h_builder( & coeff_fn, & aca );
auto A = h_builder.build( bct.get(), acc );

# Building the Right-hand Side

Building the right-hand side $$f$$ is again performed using quadrature rules over the triangular grid. The corresponding class implementing the quadrature formula is TQuadBEMRHS.

TMyRHS rhsfn;
rhs_build( 4 );
auto rhs = rhs_build.build( test_fnspace.get(), & rhsfn );

The function $$\mathbf{f}$$ is hereby provided in the form of a TBEMFunction, or, to be precise a derived class where the method eval has to be overloaded:

class TMyRHS : public TBEMFunction< complex >
{
public:
virtual complex
eval ( const T3Point & pos,
const T3Point & ) const
{
return sin( pos.x() ) * cos( pos.y() );
}
};
void eval(const TMatrix *A, TVector *v, const matop_t op, const eval_option_t &options)
evaluate LยทUยทx=y with lower triangular L and upper triangular U

In both cases, the quadrature formula and the BEM function, the value type complex and the function space (for TQuadBEMRHS) are defined as template arguments.

To bring the RHS into the H-ordering, the vector has to be permuted:

ct2->perm_e2i()->permute( rhs.get() );

# Solving the Discretised System

As standard iteration schemes will usually fail with the above equation system, H-LU preconditioning is used to ensure and to speed up convergence.

auto B = A->copy();
auto A_inv = factorise_inv( B.get(), acc );
std::unique_ptr< TFacInvMatrix > factorise_inv(TMatrix *A, const TTruncAcc &acc, const fac_options_t &options=fac_options_t())
compute factorisation of A and return inverse operator

Since the matrix is modified during LU factorisation, a copy of it has to be created and provided for factorisation. The result of factorise_inv is a linear operator suitable for evaluation of the inverse of $$A$$ and can be used for preconditioning:

TSolverInfo solve_info;
auto x = A->col_vector();
solve( A.get(), x.get(), rhs.get(), A_inv.get(), & solve_info );
std::enable_if_t< is_matrix< T1 >::value &&is_vector< T2 >::value &&is_same_type< typename T1::value_t, typename T2::value_t >::value, void > solve(T1 &A, T2 &b)
solve Aยทx = b with known A and b; x overwrites b (A is overwritten upon exit!)
Definition: Algebra.hh:990

Upon exit, x contains the computed solution to the initial discrete problem.

The ordering of the unknowns in the solution vector follows the H-ordering. To bring it back into the original ordering in the grid, use:

ct1->perm_i2e()->permute( x.get() );

The standard finalisation and catch block finishes the example:

DONE();
}
catch ( Error & e )
{
cout << e.to_string() << endl;
}
return 0;
}

# The Plain Program

#include <iostream>
#include <hlib.hh>
using namespace HLIB;
using namespace std;
class TMyRHS : public TBEMFunction< complex >
{
public:
virtual complex
eval ( const T3Point & pos,
const T3Point & ) const
{
return sin( pos.x() ) * cos( pos.y() );
}
};
int main ( int argc, char ** argv )
{
try
{
INIT();
auto grid = read_grid( "grid.tri" );
auto ansatz_fnspace = make_unique< TLinearFnSpace >( grid.get() );
auto test_fnspace = make_unique< TLinearFnSpace >( grid.get() );
auto coord1 = ansatz_fnspace->build_coord();
auto coord2 = test_fnspace->build_coord();
TAutoBSPPartStrat part_strat;
TBSPCTBuilder ct_builder( & part_strat );
auto ct1 = ct_builder.build( coord1.get() );
auto ct2 = ct_builder.build( coord2.get() );
complex kappa( -2, 1 );
TBCBuilder bct_builder;
auto bct = bct_builder.build( ct1.get(), ct2.get(), & adm_cond );
using slpbf_t = THelmholtzSLPBF< TLinearFnSpace, TLinearFnSpace >;
slpbf_t slpbf( kappa, ansatz_fnspace.get(), test_fnspace.get() );
TBFCoeffFn< slpbf_t > slp_coeff_fn( & slpbf );
TPermCoeffnFn< complex > coeff_fn( & slp_coeff_fn, bct->row_ct()->perm_i2e(), bct->col_ct()->perm_i2e() );
TACAPlus< complex > aca( & coeff_fn );
TTruncAcc acc( 1e-4, 0.0 );
TDenseMBuilder< complex > h_builder( & coeff_fn, & aca );
auto A = h_builder.build( bct.get(), acc );
TMyRHS rhsfn;
rhs_build( 4 );
auto rhs = rhs_build.build( test_fnspace.get(), & rhsfn );
ct2->perm_e2i()->permute( rhs.get() );
auto B = A->copy();
auto A_inv = factorise_inv( B.get(), acc );
TSolverInfo solve_info;
auto x = A->col_vector();
solve( A.get(), x.get(), rhs.get(), A_inv.get(), & solve_info );
ct1->perm_i2e()->permute( x.get() );
DONE();
}
catch ( Error & e )
{
cout << e.to_string() << endl;
}
return 0;
}
std::unique_ptr< TGrid > read_grid(const std::string &filename)
Read grid from file with automatic file format detection.