Python API

Python API#

The BEM interface follows the usual NGSolve form language and mirrors the mathematical notation for boundary integral operators.

Imports#

from ngsolve import *
from ngsolve.bem import *

Boundary Spaces#

Common scalar trace spaces:

  • SurfaceL2: boundary density space, used for Neumann type traces,

  • H1: volume space whose boundary trace is used for Dirichlet type traces,

  • vector valued operators use spaces such as VectorH1, HDivSurface, and HCurl traces,

  • triangular and quadrilateral surface elements are supported.

Kernels and Layer Potentials#

Laplace kernel:

\[ G_0(x,y) = \frac{1}{4\pi |x-y|}. \]

Helmholtz outgoing kernel:

\[ G_\kappa(x,y) = \frac{e^{i\kappa |x-y|}}{4\pi |x-y|}. \]

Single layer potential:

\[ (V\rho)(x) = \int_\Gamma G(x,y)\,\rho(y)\,d\sigma_y. \]

Double layer potential:

\[ (K\mu)(x) = \int_\Gamma \frac{\partial G(x,y)}{\partial n_y}\,\mu(y)\,d\sigma_y. \]

This is the naming convention behind LaplaceSL, LaplaceDL, HelmholtzSL, and HelmholtzDL.

From Paper to Code#

A single layer bilinear form

\[ \langle V\rho_h, \eta_h \rangle_\Gamma = \int_\Gamma \underbrace{\left( \int_\Gamma \,G_\kappa(x,y)\,\varphi_j(y) \,d\sigma_y \right) }_{\displaystyle{LaplaceSL(u*ds)}} \eta_i(x) \,d\sigma_x. \]

is written as

V = LaplaceSL(u * ds) * v * ds

u * ds marks the source density on the boundary, LaplaceSL(...) applies the layer potential, and multiplication by v * ds tests the result on the boundary.

Laplace Operators#

Single layer operator:

V = LaplaceSL(u * ds) * v * ds

Double layer operator:

K = LaplaceDL(u_h1 * ds) * v * ds

Helmholtz Operators#

kappa = 4.0

V = HelmholtzSL(u * ds, kappa) * v * ds
K = HelmholtzDL(u * ds, kappa) * v * ds
C = HelmholtzCF(u * ds, kappa) * v * ds

  • HelmholtzSL: single layer operator,

  • HelmholtzDL: double layer operator,

  • HelmholtzCF: combined field operator for formulations such as Brakhage Werner.

Lamé#

Lamé single layer potential:

\[ (V_{E,\nu}\mathbf t)(x) = \int_\Gamma \mathbf U_{E,\nu}(x,y)\,\mathbf t(y)\,d\sigma_y, \]
\[ \mathbf U_{E,\nu}(x,y) = \frac{1+\nu}{8\pi E(1-\nu)|x-y|} \left((3-4\nu)I + \frac{(x-y)(x-y)^T}{|x-y|^2}\right). \]

The matching NGSolve expression is the same single layer pattern, but with a vector valued density:

fes = VectorH1(mesh, order=order)
u, v = fes.TnT()

L = LameSL(u * ds, E=E, nu=nu) * v * ds

Maxwell#

The NGSolve syntax builds Maxwell blocks from HDivSurface and HCurl traces. Maxwell layer potentials use the outgoing Helmholtz kernel explicitly. The Maxwell single layer reads as

\[\begin{split} \begin{aligned} V(\mathbf j)(x) &= \kappa \int_\Gamma \frac{1}{4\pi}\frac{e^{i\kappa\|x-y\|}}{\|x-y\|}\, \mathbf j(y)\,d\sigma_y + \frac1\kappa \nabla_x \int_\Gamma \frac{1}{4\pi}\frac{e^{i\kappa\|x-y\|}}{\|x-y\|}\, \operatorname{div}_\Gamma \mathbf j(y)\,d\sigma_y,\\ \langle Vj,\rho \rangle &= \kappa \int_\Gamma \int_\Gamma \frac{1}{4\pi}\frac{e^{i\kappa\|x-y\|}}{\|x-y\|}\, \mathbf j(y)\,d\sigma_y \ d\sigma_x - \frac1\kappa \int_\Gamma \int_\Gamma \frac{1}{4\pi}\frac{e^{i\kappa\|x-y\|}}{\|x-y\|}\, \operatorname{div}_\Gamma \mathbf j(y)\, \operatorname{div}_\Gamma \mathbf \rho(x)\,d\sigma_y \,d\sigma_x,\\ \end{aligned} \end{split}\]

Single layer potential:

SL = kappa * HelmholtzSL(uHDiv.Trace()*ds, kappa) + 1/ kappa * grad(HelmholtzSL(div(uHDiv.Trace())*ds, kappa))

Single layer operator matrix:

V1 = HelmholtzSL( uHDiv.Trace()*ds, kappa) * vHDiv.Trace() * ds
V2 = HelmholtzSL( div(uHDiv.Trace()) * ds, kappa) * div(vHDiv.Trace()) * ds
V = kappa * V1.mat - 1/kappa * V2.mat

The double layer potential and weak form reads as

\[\begin{split} \begin{aligned} K(n\times\mathbf m)(x) &= \nabla_x\times\int_\Gamma \frac{1}{4\pi}\frac{e^{i\kappa\|x-y\|}}{\|x-y\|}\, \big(n(y)\times\mathbf m(y)\big)\,d\sigma_y, \\ \langle K(n\times\mathbf m)(x), \psi \rangle &= \int_\Gamma \int_\Gamma \nabla_x\times \frac{1}{4\pi}\frac{e^{i\kappa\|x-y\|}}{\|x-y\|}\, \big(n(y)\times\mathbf m(y)\big)\, \mathbf \psi(x)\, d\sigma_y \,d\sigma_x. \end{aligned} \end{split}\]

Double layer potential:

DL = MaxwellDL(Cross(n, m.Trace()) * ds, kappa)
# or
DL = curl(HelmholtzSL(Cross(n,m.Trace()) * ds, kappa))

Double layer operator:

K = MaxwellDL(Cross(n, m.Trace()) * ds, kappa) * w.Trace() * ds

Available building blocks: LaplaceSL, LaplaceDL, HelmholtzSL, HelmholtzDL, HelmholtzCF, LameSL, and MaxwellDL.

Potential Evaluation#

A potential operator can be evaluated from a grid function density:

SL = LaplaceSL(u * ds)
potential = SL(gfu)

Evaluation on a target boundary region can use a local expansion (FMM evaluation):

potential_on_screen = SL(gfu, target_boundary)

Differentiated potentials use the operator interface:

grad_potential = grad(SL)(gfu)