An indirect collocational Helmholtz BEM example

Table of Contents

• Introduction
• Theory
• Program structure
• The C++ code
• The Matlab code

Introduction

This tutorial demonstrates the usage of NiHu for solving an exterior Dirichlet problem in 3D by means of an indirect collocational type boundary element method.

Theory

The indirect BEM, applied to an exterior Dirichlet problem, expresses the radiated pressure from a vibrating surface by a double layer potential:

\( \displaystyle p({\bf x}) = \left(\mathcal{M}\sigma\right)_S({\bf x}) \quad {\bf x} \in F \)

where \( \sigma({\bf y}) \) denotes an appropriate surface source density function without direct physical meaning.

If the point \( \bf x \) approaches the boundary surface \( S \), a boundary integral equation is obtained:

\( \displaystyle p({\bf x}) = \left(\left[\mathcal{M} + \frac{1}{2}\mathcal{I}\right]\sigma\right)_S({\bf x}) \quad {\bf x} \in S \)

The boundary integral equation can be used to obtain the source density function on the surface, once a prescribed surface pressure is defined. In a second step, the radiated pressure is computed by evaluating the double layer potential integrals for external field points.

Program structure

We are going to implement a Matlab-C++ NiHu application, where

• the C++ executable is responsible for assembling the system matrices from the boundary surface mesh and the field point mesh
• the Matlab part defines the meshes, calls the C++ executable, solves the systems of equations, computes the radiated pressure and quantifies the error of the solution.

The C++ code

The C++ code is going to be called from Matlab as

[Ms, Mf] = helmholtz_ibem_3d(
    surf_nodes, surf_elements, field_nodes, field_elements, wave_number);

Note

The computed system matrix Ms will contain the discretised identity operator too.

Without going into details regarding the C++ code, we examplify the lines where the operators are discretised using collocation:

    double k = *mxGetPr(rhs[4]);
    auto I = NiHu::identity_integral_operator();
    auto M = NiHu::create_integral_operator(NiHu::helmholtz_3d_DLP_kernel<double>(k));
 
    Ms << dirac(surf_sp) * M[surf_sp] + dirac(surf_sp) * (.5*I)[surf_sp];
    Mf << field_sp * M[surf_sp];

The Matlab code

The example creates a sphere surface of radius \( R = 1\,\mathrm{m} \) centred at the origin, consisting of triangular elements only. The field point mesh is a surface obtained by revolving a line around the radiator.

The Dirichlet boundary conditions are defined by a point source located at the point \( \mathbf{x}_0 \).

The compiled mex file can be called from Matlab, as demonstrated in the example below:

%// sphere radiator with radius R = 1, divided into 6 elements along the radius.
radiator = quad2tria(create_sphere_boundary(1, 5));
%// field point mesh, a line revolved around the z axis
field = revolve_mesh(create_line([1.125, 0, 0; 3.625, 0, 0], 40), pi/100, 50, [0 0 1]);
 
%// extract the mesh description matrices
[r_nodes, r_elem] = extract_core_mesh(radiator);
[f_nodes, f_elem] = extract_core_mesh(field);
 
%// call C++ code at wave number k = 0
k = 4;
[Ms, Mf] = helmholtz_bem_indirect_dirichlet(r_nodes, r_elem, f_nodes, f_elem, k);
 
%// define the incident velocity field and analytical solutions
x0 = [.4 .3 0];
[r_cent, r_norm] = centnorm(radiator);
[ps_ana, qs_ana] = incident('point', x0, r_cent, r_norm, k);
 
%// analytical solution on the field
f_cent = centnorm(field);
pf_ana = incident('point', x0, f_cent, [], k);
 
%// acoustic pressure on the surface and in the field points
sigma_inf = Ms \ ps_ana;
pf_num = Mf * sigma_inf;
 
% // calculate errors
pf_err = abs(pf_num./pf_ana - 1);
 
% // plot results
figure;
subplot(1,2,1);
plot_mesh(field, real(pf_num)); view(2);  shading flat;
c1 = colorbar('SouthOutside');
xlabel(c1, 'Real part of numerical solution');
subplot(1,2,2);
plot_mesh(field, log10(pf_err)); view(2);  shading flat;
c1 = colorbar('SouthOutside');
xlabel(c1, 'Log10 error');
 
%// display error information
fprintf('Field   log10 error: %f \n', log10(mean(pf_err)));

The generated plot looks like

And the resulting errors read as

log10 mean error on field = -2.42

The full codes of this tutorial are available here:

• C++ code: tutorial/helmholtz_ibem_3d.mex.cpp
• Matlab code: tutorial/helmholtz_ibem_3d_test.m