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 );
THiLoFreqGeomAdmCond       adm_cond( kappa, 10 );
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;
TQuadBEMRHS< TLinearFnSpace, complex >  
                           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() );
    }
};

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 );

HLIB::factorise_inv

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 );

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 );
    THiLoFreqGeomAdmCond       adm_cond( kappa, 10 );
    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;
    TQuadBEMRHS< TLinearFnSpace, complex >
                               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;
}