Helmholtz#

Goal#

Solve a Helmholtz boundary integral equation (frequency-domain acoustics) on a surface mesh.

This tutorial corresponds to the example program:

  • examples/helmholtz_bem.cpp

What the code does#

  1. Loads a surface mesh.

  2. Sets the wavenumber k.

  3. Builds the Helmholtz boundary integral operator(s).

  4. Assembles the right-hand side.

  5. Solves with GMRES and writes VTK output.

How to run#

./build/helmholtz_bem ../sphere_0.1.off 5.0

Full source#

 1#include "biekit/biekit.h"
 2
 3#include <complex>
 4#include <iostream>
 5#include <string>
 6#include <vector>
 7
 8int main(int argc, char** argv) {
 9    std::string mesh_path = "../sphere_0.1.off";
10    if (argc >= 2) {
11        mesh_path = argv[1];
12    }
13
14    std::cout << "Loading mesh: " << mesh_path << "\n";
15    Mesh mesh = load_mesh(mesh_path);
16
17    const std::array<double, 3> khat{0.0, 0.0, 1.0};
18    const double wavelength = 1.0;
19    const std::complex<double> k = std::complex<double>{2.0 * kPi / wavelength, 0.0};
20
21    ScalarBemOptions opt;
22    opt.quad_order = 2;
23    opt.enable_singularity_correction = true;
24
25    ScalarSingleLayerOperator A(mesh, k, opt);
26
27    std::vector<std::complex<double>> rhs = assemble_dirichlet_rhs_plane_wave(mesh, k, khat);
28    std::vector<std::complex<double>> sol(rhs.size(), std::complex<double>{0.0, 0.0});
29
30    std::cout << "Solving Helmholtz BEM with GMRES...\n";
31    GmresComplex solver;
32    const auto info = solver.solve(A, sol, rhs);
33    std::cout << "GMRES converged=" << (info.converged ? "true" : "false")
34              << "  rel_res=" << info.rel_residual << "\n";
35
36    std::vector<double> sigma_mag(sol.size(), 0.0);
37    for (std::size_t i = 0; i < sol.size(); ++i) {
38        sigma_mag[i] = std::abs(sol[i]);
39    }
40
41    VtkWriteOptions vtk_opt;
42    vtk_opt.title = "helmholtz_bem_dirichlet";
43    WriteVtkPolyDataCellFields(
44        mesh,
45        "helmholtz_solution.vtk",
46        {"sigma_mag"},
47        {sigma_mag},
48        vtk_opt);
49
50    std::cout << "[success] write vtk: helmholtz_solution.vtk\n";
51    return info.converged ? 0 : 2;
52}