Laplace#

Goal#

Solve a Laplace boundary integral equation on a closed surface mesh.

This tutorial corresponds to the example program:

  • examples/laplace_bem.cpp

What the code does#

  1. Loads a surface mesh.

  2. Builds the Laplace operators (single-/double-layer, depending on the example implementation).

  3. Assembles the right-hand side from a simple manufactured solution / point source.

  4. Solves the linear system (GMRES).

  5. Writes VTK output for visualization.

How to run#

./build/laplace_bem ../sphere_0.1.off

Full source#

 1#include "biekit/biekit.h"
 2
 3#include <cmath>
 4#include <complex>
 5#include <iostream>
 6#include <string>
 7#include <vector>
 8
 9int main(int argc, char** argv) {
10    std::string mesh_path = "../sphere_0.1.off";
11    if (argc >= 2) {
12        mesh_path = argv[1];
13    }
14
15    std::cout << "Loading mesh: " << mesh_path << "\n";
16    Mesh mesh = load_mesh(mesh_path);
17
18    const std::complex<double> k_laplace{0.0, 0.0};
19
20    ScalarBemOptions opt;
21    opt.quad_order = 2;
22    opt.enable_singularity_correction = true;
23
24    ScalarSingleLayerOperator A(mesh, k_laplace, opt);
25
26    const std::array<double, 3> src_point{0.0, 0.0, 2.0};
27    const auto centers = face_centers_tri(mesh);
28    std::vector<std::complex<double>> rhs;
29    rhs.resize(centers.size());
30
31    for (std::size_t i = 0; i < centers.size(); ++i) {
32        const auto& r = centers[i];
33        const double dx = r[0] - src_point[0];
34        const double dy = r[1] - src_point[1];
35        const double dz = r[2] - src_point[2];
36        const double dist = std::sqrt(dx * dx + dy * dy + dz * dz);
37        rhs[i] = -std::complex<double>{LaplaceGreenFunction(dist), 0.0};
38    }
39    std::vector<std::complex<double>> sol(rhs.size(), std::complex<double>{0.0, 0.0});
40
41    std::cout << "Solving Laplace BEM with GMRES...\n";
42    GmresComplex solver;
43    const auto info = solver.solve(A, sol, rhs);
44    std::cout << "GMRES converged=" << (info.converged ? "true" : "false")
45              << "  rel_res=" << info.rel_residual << "\n";
46
47    std::vector<double> sigma_mag(sol.size(), 0.0);
48    for (std::size_t i = 0; i < sol.size(); ++i) {
49        sigma_mag[i] = std::abs(sol[i]);
50    }
51
52    VtkWriteOptions vtk_opt;
53    vtk_opt.title = "laplace_bem_dirichlet";
54    WriteVtkPolyDataCellFields(
55        mesh,
56        "laplace_solution.vtk",
57        {"sigma_mag"},
58        {sigma_mag},
59        vtk_opt);
60
61    std::cout << "[success] write vtk: laplace_solution.vtk\n";
62    return info.converged ? 0 : 2;
63}