Feature Request: Stabilized Maxwell BEM for Low-Frequency Problems

[ksugahar]

Feature Request: Stabilized Maxwell BEM for Low-Frequency Problems

Background

This issue follows a discussion in the ngbem repository (Weggler/ngbem#14) with @Weggler (Lucy Weggler), who suggested implementing the stabilized Maxwell formulation in NGSolve core rather than the legacy ngbem addon.

The stabilized formulation is based on:

L. Weggler, "High Order Boundary Element Methods,"
Ph.D. Dissertation, Universitaet des Saarlandes, 2011.

• Chapter 3.3: Stabilized Formulation (Theory)
• Chapter 5.3: Stabilized Formulation (Discretization)

Problem: Low-Frequency Breakdown

The classical EFIE (Electric Field Integral Equation) formulation suffers from low-frequency breakdown:

cond(A_c) ~ O(kappa^{-2})  as kappa -> 0

From the thesis (Table 6.6), for a sphere with N=128 elements, p=1:

kappa [1/m]	Classical cond(A_c)	Stabilized cond(A_s)
5.0	3.9e+1	1.6e+1
5.0e-1	3.5e+2	1.1e+2
5.0e-2	3.5e+4	6.2e+1
5.0e-3	3.5e+6	6.2e+1
5.0e-4	3.5e+8	6.2e+1

This limits the applicability of Maxwell BEM to high-frequency problems only.

Solution: Stabilized Formulation

The stabilized formulation (Eq. 5.20) introduces an additional unknown (surface charge density rho_Gamma) and explicitly recovers the Gauss law:

| A_kappa    Q_kappa      |   | j^t       |   | M*m |
|                         | * |           | = |     |
| Q_kappa^T  kappa^2*V    |   | rho_Gamma |   |  0  |

where:

• A_kappa: Maxwell single layer potential (already in NGSolve)
• V_kappa: Helmholtz single layer potential (already in NGSolve)
• Q_kappa: Transition matrix (new - Eq. 5.19)

The transition matrix couples the surface divergence of HDivSurface elements with the scalar H1 space:

(Q_kappa)_lk = int_Gamma int_Gamma div_Gamma phi_l(x) G_kappa(x-y) nu_k(y) dsigma_y dsigma_x

What is Currently Available

In ngsolve/bem/:

• MaxwellSLKernel<3>, MaxwellDLKernel<3> in kernels.hpp
• DiffOpMaxwell, DiffOpMaxwellNew in bem_diffops.hpp
• Python bindings: MaxwellSingleLayerPotentialOperator, MaxwellDoubleLayerPotentialOperator

What is Missing

As @Weggler noted, the missing piece is:

"there has been only one operator missing, right? The one that takes the normal trace of E (the normal trace of the potential)"

Specifically, we need:

1. A way to construct the transition operator Q_kappa
2. This requires coupling div_Gamma (from HDivSurface) with the scalar Helmholtz single layer potential

Proposed Implementation

Option 1: New Kernel + DiffOp

Add a TransitionKernel<3> that uses the Helmholtz Green's function with a trial evaluator computing div_Gamma:

template<>
class TransitionKernel<3> : public BaseKernel
{
    double kappa;
public:
    typedef Complex value_type;
    static string Name() { return "Transition"; }

    template <typename T>
    auto Evaluate(Vec<3,T> x, Vec<3,T> y, Vec<3,T> nx, Vec<3,T> ny) const
    {
        T r = L2Norm(x - y);
        auto G = exp(Complex(0, kappa) * r) / (4.0 * M_PI * r);
        return Vec<1, Complex>(G);
    }
    // ... FMM support similar to HelmholtzSLKernel
};

Option 2: Compose Existing Operators

Use the relationship from Eq. 5.22:

Q_kappa = D * V_kappa

where D is the discrete surface divergence operator.

References

• Original discussion: Feature Proposal: Stabilized Maxwell BEM for Low-Frequency Problems (based on your thesis) Weggler/ngbem#14
• Prototype code (Gist): https://gist.github.com/ksugahar/985063de1b5a1abea040f2c744b7d61b
• Thesis reference:
  • Eq. 3.77: Continuity equation (weak form)
  • Eq. 5.14: Stabilized variational formulation
  • Eq. 5.19: Transition matrix definition
  • Eq. 5.20: Block system matrix
  • Lemma 5.3.1: Equivalence proof via Schur complement

Deliverables

If this feature is accepted, I am happy to contribute:

1. Implementation of the transition operator
2. Python bindings
3. Documentation notebook demonstrating the stabilized formulation
4. Test cases comparing classical vs stabilized at various kappa values

Looking forward to your feedback!

cc: @Weggler @JSchoeberl

[ksugahar]

Additional Approach: Loop-Star Decomposition

In addition to the stabilized formulation discussed above, there is another well-established approach for low-frequency EFIE stabilization: Loop-Star decomposition.

Principle

The Loop-Star decomposition separates the current into:

• Loop basis: Solenoidal component (kernel of div operator)
• Star basis: Non-solenoidal component (range of grad operator)

In this basis, the EFIE matrix becomes block-diagonal at DC:

Z_LS = | Z_LL   0    |
       | 0     Z_SS  |

The condition number issue is isolated to the Star-Star block, which can be handled by appropriate scaling (1/omega^2).

References

1. Vecchi, G. "Loop-star decomposition of basis functions in the discretization of the EFIE", IEEE Trans. Antennas Propagat., Vol. 47, No. 2, pp. 339-346, 1999.

2. Meddahi, S. & Selgas, V. "A mixed-FEM and BEM coupling for the approximation of the scattering of thermal waves in locally non-homogeneous media", ESAIM: M2AN, Vol. 35, No. 6, 2001.

Comparison with Stabilized Formulation

Approach	Additional Unknowns	Basis Transformation	High-Order Extension
Stabilized (Weggler)	Surface charge rho_Gamma	No	Natural
Loop-Star	None	Yes (spanning tree)	Requires careful design

Question

Is there interest in supporting both approaches in NGSolve?

The stabilized formulation (discussed in the original issue) is mathematically elegant and extends naturally to high-order elements. The Loop-Star approach is widely used in CEM community and may be more familiar to some users.

I have a prototype Loop-Star implementation available if there is interest.