Table of contents

In this example, the scattering of an incoming plane wave between two bodies is simulated as depicted in the following figure:

Scattering of an incoming plane wave

The mathematical formulation leads to an integral equation based on the Helmholtz kernel. Solving such an equation involves the discrete representation of the analytic operators and, if neccessary, the computation of a preconditioner. Both shall be performed using the 𝓗-matrix technique.

Problem Description

Let $\Omega_i \in R^3, i = 1, \ldots, r$ be convex sets with smooth boundaries $\Gamma_i = \partial \Omega_i$ . Furthermore, let $\nu = \nu_1, \nu_2, \ldots$ with $1 \le \nu_i \le r$ be a sequence, such that $\nu_i \ne \nu_{i+1}$ for all $i$ .

To compute: densities $\eta_{\nu_i}$ on $\Gamma_{\nu_i}$ defined by

$\eta_{\nu_i} - \mathcal{K}_{\nu_i} \eta_{\nu_i} = g_{\nu_i},$

with

$(\mathcal{K}_{\nu_i} \eta_{\nu_i})(x) := \int_{\Gamma_{\nu_i}} \frac{\partial G(x,y)}{\partial n(x)} \eta_{\nu_i}(y) dy, \quad x \in \Gamma_{\nu_i}$

with $G(x,y) = -\frac{1}{2\pi} \frac{e^{ik|x-y|}}{|x-y|}$ , $k > 0$ and with $d \in R^3,\; |d| = 1,$

$\begin{eqnarray*} g_{\nu_i}(x) := \left\{ \begin{array}{ll} 2 \frac{\partial e^{ik \langle d,x \rangle }}{\partial n(x)}, & i = 1 \\ (\mathcal{K}_{\nu_{i-1}} \eta_{\nu_{i-1}})(x), \; x \in \Gamma_{\nu_i} & i > 1 \end{array} \right. \end{eqnarray*}$

With a standard Galerkin discretisation using constant ansatz functions, one obtains the following matrices

$K^{\nu_i,\nu_i}_{ij} := \frac{1}{2\pi} \int_{\Gamma_{\nu_i}} \int_{\Gamma_{\nu_i}} e^{ik| x - y|} \frac{1-ik|x-y|}{|x-y|^3} \langle x-y, n( x ) \rangle dy dx$

and

$K^{{\nu_i},{\nu_{i-1}}}_{ij} := \frac{1}{2\pi} \int_{\Gamma_{\nu_i}} \int_{\Gamma_{\nu_{i-1}}} e^{ik| x - y|} \frac{1-ik|x-y|}{|x-y|^3} \langle x-y, n( x ) \rangle dy dx$

This leads to the following algorithm, which iterates between the surfaces $v_i$ :

build $K^{1,1}, K^{2,2}, K^{2,1}$ and $K^{1,2}$ ;
build initial right hand side $b_{1,1} := g_1$ ;
for
- solve $(I - K^{1,1}) \eta_{i,1} = b_{i,1}$ ;
- if $\eta_{i,1} / \eta_{i-1,1}$ converged stop
- $b_{i,2} := K^{2,1} \eta_{i,1}$ ;
- solve $(I - K^{2,2}) \eta_{i,2} = b_{i,2}$ ;
- $b_{i+1,1} := K^{1,2} \eta_{i,2}$ ;

Grid Definition

After the standard initialisation, which also imports the 𝓗𝖫𝖨𝖡𝗉𝗋𝗈 types for real and complex values,

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

the program begins with the definition of the grid for the surfaces $\Gamma_1$ and $\Gamma_2$ . While the grid is read from a file, we assume, that it contains a sphere, centered at $(0,0,0)$ with radius 1.

auto                    grid1( read_grid( "sphere.tri" ) );
grid1->scale( T3Point( 4, 4, 1 ) );

The second command is optional. It deforms the sphere into a disc, which extends the scattering process.

Usually, a grid in 𝓗𝖫𝖨𝖡𝗉𝗋𝗈 will only define triangle normals. Since this leads to constant normal directions per triangle, the approximation of the Helmholtz DLP kernel is less accurate and may even lead to false results (think of scattering between parallel planes). Therefore, the grid file should contain vertex normals, which are used to interpolate normal directions for each point on the triangle, leading to a continuous normal direction on the surface.

If vertex normals are missing, the can be trivially added under the above assumptions since each vertex coordinate equals the corresponding normal direction. You may use add_vtx_normal of TGrid for this.

Since we need two surfaces, the process is repeated, but in addition, the second disc is shifted with respect to the z direction, such that both disc will have a positive distance:

auto                    grid2( read_grid( "sphere.tri" ) );
grid2->scale( T3Point( 4, 4, 1 ) );
grid2->translate( T3Point( 0, 0, 2.5 ) );  // distance 0.5

This leads to the following configuration:

Grid with two discs

The construction of the grid ends with the definition of the corresponding function spaces. Here, we restrict ourselves to piecewise constant ansatz:

TConstFnSpace fnspace1( grid1.get() );

TConstFnSpace fnspace2( grid2.get() );

Matrix Construction

We need to construct the matrices $K^{1,1}, K^{2,2}, K^{2,1}$ and $K^{1,2}$ . Due to definition of the grid, the matrices $K^{1,1}$ and $K^{2,2}$ are identical.

For the remaining matrices, cluster trees and block cluster trees are needed:

auto                    coord1 = fnspace1.build_coord();
auto                    coord2 = fnspace2.build_coord();
TAutoBSPPartStrat       part_strat;
TBSPCTBuilder           ct_builder( & part_strat );
auto                    ct1 = ct_builder.build( coord1.get() );
auto                    ct2 = ct_builder.build( coord2.get() );
TStdGeomAdmCond         adm_cond( 2.0 );
TBCBuilder              bct_builder;
auto                    bct11 = bct_builder.build( ct1.get(), ct1.get(), & adm_cond );
auto                    bct12 = bct_builder.build( ct1.get(), ct2.get(), & adm_cond );
auto                    bct21 = bct_builder.build( ct2.get(), ct1.get(), & adm_cond );

First, the matrix $K^{1,1}$ shall be constructed. 𝓗𝖫𝖨𝖡𝗉𝗋𝗈 provides the class TAcousticScatterBF, which implements the kernel for $G(x,y)$ . Not defined so far is the wavenumber $\kappa$ , which should have a value of 6. Furthermore, ACA+ is used to compute low-rank approximations with a fixed block-wise accuracy of $10^{-4}$ .

using  bf = TAcousticScatterBF< TConstFnSpace, TConstFnSpace >;
complex_t                    kappa( 4, 0 );
TTruncAcc                    acc = fixed_prec( 1e-4 );
unique_ptr< TMatrix >        K11;
{
    bf_t                         bf( kappa, & fnspace1, & fnspace1 );
    TBFCoeffFn< bf_t >           coeff_fn( & bf, bct11->row_ct()->perm_i2e(), bct11->col_ct()->perm_i2e() );
    TACAPlus< complex_t >        lrapx( & coeff_fn );
    TDenseMBuilder< complex_t >  h_builder( & coeff_fn, & lrapx );
        
    K11 = unique_ptr< TMatrix >( h_builder.build( bct11.get(), acc ) );
}

Remarks: Objects such as the bilinear form, the coefficient function or the low-rank approximation, which are specific for a single matrix are encapsulated in a block to prevent name clashes with the other matrices.

The actual operator in the above equation is $I - K^{1,1}$ . Here $I$ is not the identity, but the mass matrix, as defined by the identity kernel:

using  massbf_t = TMassBF< TConstFnSpace, TConstFnSpace >;
{
    massbf_t                     bf( & fnspace1, & fnspace1 );
    TBFCoeffFn< massbf_t >       coeff_fn( & bf, bct11->row_ct()->perm_i2e(), bct11->col_ct()->perm_i2e() );
    TACAPlus< real_t >           lrapx( & coeff_fn );
    TDenseMBuilder< real_t >     h_builder( & coeff_fn, & lrapx );
    
    M11 = unique_ptr< TMatrix >( h_builder.build( bct11.get(), acc ) );
    M11->set_complex( true );
}

Finally, $K^{1,1}$ is updated with the mass matrix:

add( 1.0, M11.get(), -1.0, K11.get(), acc );

The construction of $K^{1,2}$ and $K^{2,1}$ is similar to $K^{1,1}$ . Only the function spaces and block cluster trees will change.


$K^{1,1}$	$K^{1,2}$	$K^{2,1}$

Please note, that $K^{1,1}$ is more dense than the other two matrices. This is due to the positive distance between both discs. If this is increased (or if coarsening is applied), $K^{1,2}$ and $K^{2,1}$ may be represented by a single low-rank matrix.

Right-Hand Side

The initial right-hand side is defined by $g_{\nu_i}$ above, where also the the incident field $d$ is used:

$g_{\nu_i}(x) = 2 \frac{\partial e^{ik \langle d,x \rangle }}{\partial n(x)}$

Under the above setting, this is equal to:

$g_{\nu_i}(x) = 2 i \kappa e^{ ik \langle d,x \rangle } \langle d, n \rangle$

and is implemented in the following RHS function:

class TAcousticRHS : public TBEMFunction< complex_t >
{
private:
    const T3Point    _inc_dir;  // incident direction
    const complex_t  _ikappa;
    
public:
    TAcousticRHS ( const T3Point &  inc_dir,
                   const complex_t  kappa )
            : _inc_dir( inc_dir )
            , _ikappa( complex_t( 0, 1 ) * kappa )
    {}
    
    virtual complex_t
    eval ( const T3Point &  pos,
           const T3Point &  normal ) const
    {
        return  real_t(2) * _ikappa * exp( _ikappa * dot( _inc_dir, pos ) ) * dot( _inc_dir, normal );
    }
};

We use TQuadBEMRHS to compute the corresponding vector of the RHS using quadrature and afterwards reorder it for 𝓗-arithmetic:

TAcousticRHS                     rhsfn( T3Point( 1, 1, 0 ), kappa );
TQuadBEMRHS< TConstFnSpace, complex_t >  
                                 rhs_build( 4 );
auto                             rhs = rhs_build.build( & fnspace1, & rhsfn );
ct1->perm_e2i()->permute( rhs.get() );

The plane wave $d$ is chosen to be $(1,1,0)$ .

Acoustic Scattering Iteration

The scattering iteration involves to solves with $K^{1,1}$ ( $=K^{2,2}$ ). This is accelerated using a preconditioner computed by and 𝓗-LU factorisation:

auto K11_copy = K11->copy();

auto K11_pre = factorise_inv( K11_copy.get(), acc );

We also need a vector to hold the current and the last approximation of the solution:

auto                             x     = K11->col_vector();
auto                             x_old = K11->col_vector();
real_t                           old_ratio = real_t(1);

The variable old_ratio is needed to check convergence.

The iteration begins with solving $(I - K^{1,1}) \eta_{i,1} = b_{i,1}$ , where $b$ is the vector of the RHS:

for ( int it = 0; it < 100; ++it )
{
    solve( K11.get(), x.get(), rhs.get(), K11_pre.get() );

Convergence is tested, if the ratio between the previous solution and the current solution is below some threshold. For this ratio is tested entry-wise, e.g. $\eta_{i,1}^j / \eta_{i-1,1}^j$ for each $j$ . This is implemented in the function compare:

real_t
compare ( const TVector * x,
          const TVector * x_old )
{
    real_t  vmax = 0.0;
    size_t  n   = x->size();
    for ( idx_t  j = 0; j < idx_t(n); ++j )
    {
        const complex_t  v     = x->centry( j );
        const complex_t  v_old = x_old->centry( j );
        
        vmax = std::max( vmax, abs( v / v_old ) );
    }
    return vmax;
}

If this ratio does no longer change, the iteration is stopped. The threshold in this example is $10^{-4}$ :

if ( it > 0 )
{
    const real_t  ratio = compare( x.get(), x_old.get() );
    const real_t  diff  = ( it > 1 ? std::abs( ratio - old_ratio ) / std::abs( old_ratio ) : 1.0 );
    old_ratio = ratio;
    if ( diff < 1e-4 )
        break;
}

In addition to this, the iteration shall stop if the solution $\eta_{i,1}$ vanishes:

const real_t  normx = x->norm2();
if ( x->norm2() < 1e-20 )
    break;

Otherwise, the iteration continues with the scaling of the solution vector and the computation of the RHS for $K^{2,2}$ :

x->scale( real_t(1) / normx );
x_old->assign( real_t(1), x.get() );
mul_vec( real_t(-1), K12.get(), x.get(), real_t(0), rhs.get(), MATOP_NORM );

Solving with $K^{2,2}$ ( $=K^{1,1}$ ) and computing the RHS for the next iteration finished the loop and the program:

        solve( K11.get(), x.get(), rhs.get(), K11_pre.get(), & solve_info );
        mul_vec( real_t(-1), K21.get(), x.get(), real_t(0), rhs.get(), MATOP_NORM );
      }
      DONE();
  }
  catch ( Error & e )
  {
      cout << e.to_string() << endl;
  }
    
  return 0;
}

Running the simulation will produce the following output for $\eta_{i,1}$ and $\eta_{i,2}$ . Here, the entry-wise absolute values of the complex valued solution vectors are used for the visualisation (reload page to restart animation):


$\eta_{i,1}$	$\eta_{i,2}$

The Plain Program

Remarks: Within the 𝓗𝖫𝖨𝖡𝗉𝗋𝗈 distribution you will find the full program with visualisation commands.

#include <iostream>
#include <hlib.hh>
using namespace HLIB;
using namespace std;
using  real_t     = HLIB::real;
using  copmplex_t = HLIB::complex;
class TAcousticRHS : public TBEMFunction< complex_t >
{
private:
    const T3Point    _inc_dir;  // incident direction
    const complex_t  _ikappa;
    
public:
    TAcousticRHS ( const T3Point &  inc_dir,
                   const complex_t  kappa )
            : _inc_dir( inc_dir )
            , _ikappa( complex_t( 0, 1 ) * kappa )
    {}
    
    virtual complex_t
    eval ( const T3Point &  pos,
           const T3Point &  normal ) const
    {
        return  real_t(2) * _ikappa * exp( _ikappa * dot( _inc_dir, pos ) ) * dot( _inc_dir, normal );
    }
};
real_t
compare ( const TVector * x,
          const TVector * x_old )
{
    real_t  vmax = 0.0;
    size_t  n    = x->size();
    for ( idx_t  j = 0; j < idx_t(n); ++j )
    {
        const complex_t  v     = x->centry( j );
        const complex_t  v_old = x_old->centry( j );
        
        vmax = std::max( vmax, abs( v / v_old ) );
    }// for
    return vmax;
}
int main ( int argc, char ** argv ) 
{
  try
  {
    INIT();
    auto                    grid1( read_grid( "sphere.tri" ) );
    grid1->scale( T3Point( 4, 4, 1 ) );
    auto                    grid2( read_grid( "sphere.tri" ) );
    grid2->scale( T3Point( 4, 4, 1 ) );
    grid2->translate( T3Point( 0, 0, 2.5 ) );  // distance 0.5
    TConstFnSpace           fnspace1( grid1.get() );
    TConstFnSpace           fnspace2( grid2.get() );
    auto                    coord1 = fnspace1.build_coord();
    auto                    coord2 = fnspace2.build_coord();
    TAutoBSPPartStrat       part_strat;
    TBSPCTBuilder           ct_builder( & part_strat );
    auto                    ct1 = ct_builder.build( coord1.get() );
    auto                    ct2 = ct_builder.build( coord2.get() );
    
    TStdGeomAdmCond         adm_cond( 2.0 );
    TBCBuilder              bct_builder;
    auto                    bct11 = bct_builder.build( ct1.get(), ct1.get(), & adm_cond );
    auto                    bct12 = bct_builder.build( ct1.get(), ct2.get(), & adm_cond );
    auto                    bct21 = bct_builder.build( ct2.get(), ct1.get(), & adm_cond );
    using  bf_t = TAcousticScatterBF< TConstFnSpace, TConstFnSpace >;
    complex_t               kappa( 6, 0 );
    TTruncAcc               acc = fixed_prec( 1e-4 );
    unique_ptr< TMatrix >   K11;
    {
        bf_t                         bf( kappa, & fnspace1, & fnspace1 );
        TBFCoeffFn< bf_t >           coeff_fn( & bf, bct11->row_ct()->perm_i2e(), bct11->col_ct()->perm_i2e() );
        TACAPlus< complex_t >        lrapx( & coeff_fn );
        TDenseMBuilder< complex_t >  h_builder( & coeff_fn, & lrapx );
            
        K11 = h_builder.build( bct11.get(), acc );
    }
    using  massbf_t = TMassBF< TConstFnSpace, TConstFnSpace >;
    unique_ptr< TMatrix >        M11;
    {
        massbf_t                     bf( & fnspace1, & fnspace1 );
        TBFCoeffFn< massbf_t >       coeff_fn( & bf, bct11->row_ct()->perm_i2e(), bct11->col_ct()->perm_i2e() );
        TACAPlus< real_t >           lrapx( & coeff_fn );
        TDenseMBuilder< real_t >     h_builder( & coeff_fn, & lrapx );
        
        M11 = h_builder.build( bct11.get(), acc );
        M11->set_complex( true );
    }
    add( 1.0, M11.get(), -1.0, K11.get(), acc );
    unique_ptr< TMatrix >        K12;
    {
        bf_t                         bf( kappa, & fnspace1, & fnspace2 );
        TBFCoeffFn< bf_t >           coeff_fn( & bf, bct12->row_ct()->perm_i2e(), bct12->col_ct()->perm_i2e() );
        TACAPlus< complex_t >        lrapx( & coeff_fn );
        TDenseMBuilder< complex_t >  h_builder( & coeff_fn, & lrapx );
            
        K12 = h_builder.build( bct12.get(), acc );
    }
    unique_ptr< TMatrix >        K21;
    {
        bf_t                         bf( kappa, & fnspace2, & fnspace1 );
        TBFCoeffFn< bf_t >           coeff_fn( & bf, bct21->row_ct()->perm_i2e(), bct21->col_ct()->perm_i2e() );
        TACAPlus< complex_t >        lrapx( & coeff_fn );
        TDenseMBuilder< complex_t >  h_builder( & coeff_fn, & lrapx );
            
        K21 = h_builder.build( bct21.get(), acc );
    }
    TAcousticRHS                     rhsfn( T3Point( 1, 1, 0 ), kappa );
    TQuadBEMRHS< TConstFnSpace, complex_t >  
                                     rhs_build( 4 );
    auto                             rhs = rhs_build.build( & fnspace1, & rhsfn );
    ct1->perm_e2i()->permute( rhs.get() );
    auto                             K11_copy = K11->copy();
    auto                             K11_pre  = factorise_inv( K11_copy.get(), acc );
    auto                             x     = K11->col_vector();
    auto                             x_old = K11->col_vector();
    real_t                           old_ratio = real_t(1);
    for ( int it = 0; it < 100; ++it )
    {
        solve( K11.get(), x.get(), rhs.get(), K11_pre.get() );
        if ( it > 0 )
        {
            const real_t  ratio = compare( x.get(), x_old.get() );
            const real_t  diff  = ( it > 1 ? std::abs( ratio - old_ratio ) / std::abs( old_ratio ) : 1.0 );
            old_ratio = ratio;
            if ( diff < 1e-4 )
                break;
        }
        const real_t  normx = x->norm2();
        if ( x->norm2() < 1e-20 )
            break;
        x->scale( real_t(1) / normx );
        x_old->assign( real_t(1), x.get() );
        mul_vec( real_t(-1), K12.get(), x.get(), real_t(0), rhs.get(), MATOP_NORM );
        solve( K11.get(), x.get(), rhs.get(), K11_pre.get() );
        mul_vec( real_t(-1), K21.get(), x.get(), real_t(0), rhs.get(), MATOP_NORM );
      }
      DONE();
  }
  catch ( Error & e )
  {
      cout << e.to_string() << endl;
  }
    
  return 0;
}