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 <hpro.hh>
 
using namespace Hpro;
using namespace std;
 
using real_t    = double;
using complex_t = std::complex< real_t >;
 
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. The function spaces permit different value types, i.e., float and double.

auto ansatz_fnspace = make_unique< TLinearFnSpace< real_t > >( grid.get() );

auto test_fnspace = make_unique< TLinearFnSpace< real_t > >( 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_t             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< real_t >, TLinearFnSpace< real_t >, complex_t >;
 
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_t > 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_t >        aca( & coeff_fn );
TTruncAcc                    acc( 1e-4, 0.0 );
TDenseMBuilder< complex_t >  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.

auto  rhsfn     = TMyRHS();
auto  rhs_build = TQuadBEMRHS< TLinearFnSpace, TMyRHS, complex_t >( 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_t >
{
public:
    virtual complex_t
    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_t 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 );

Hpro::factorise_inv

std::unique_ptr< TFacInvMatrix< value_t > > factorise_inv(TMatrix< value_t > *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 Hpro;
using namespace std;
  
using real_t    = double;
using complex_t = std::complex< real_t >;
 
class TMyRHS : public TBEMFunction< complex_t >
{
public:
    virtual complex_t
    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< real_t > >( grid.get() );
    auto                         test_fnspace   = make_unique< TLinearFnSpace< real_t > >( 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_t                    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< real_t >, TLinearFnSpace< real_t >, complex_t >;
 
    slpbf_t                      slpbf( kappa, ansatz_fnspace.get(), test_fnspace.get() );
    TBFCoeffFn< slpbf_t >        slp_coeff_fn( & slpbf );
    TPermCoeffnFn< complex_t >   coeff_fn( & slp_coeff_fn, bct->row_ct()->perm_i2e(), bct->col_ct()->perm_i2e() );
    TACAPlus< complex_t >        aca( & coeff_fn );
    TTruncAcc                    acc( 1e-4, 0.0 );
    TDenseMBuilder< complex_t >  h_builder( & coeff_fn, & aca );
    auto                         A = h_builder.build( bct.get(), acc );
 
    auto                         rhsfn     = TMyRHS();
    auto                         rhs_build = TQuadBEMRHS< TLinearFnSpace, complex_t >( 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;
}