bem — Boundary-Element Collocation

Status: \(n \le 2\) only, and numerically identical to mom

Since 0.11.0 bem rejects \(n \ge 3\) soils with NotImplementedError (it shared the structurally incomplete complex-image kernel with cim — audit 2026-07-08, WP-E), and since 0.15.0 it computes no complex-image fit at all. In the two regimes it accepts it assembles the same reaction matrix as mom and solves it with the same constraint solver; the measured relative difference on a ring case is 4e-16. bem is therefore not an independent cross-check of mom — treating "mom and bem agree" as corroboration counts one computation twice (ADR-0002 amendment 2026-07-09). It is kept as the collocation-flavoured entry point of the family.

Physical context

The Boundary-Element Method is a sister to the Method of Moments: both reduce a continuous boundary integral equation to a finite linear system by discretising the boundary into elements (here: wire segments) and choosing a finite set of test functions. The two methods only differ in the choice of test function:

• Galerkin (mom): test function = basis function. The reaction matrix entry \(Z_{ij}\) is the average potential of segment \(i\) caused by a unit current on segment \(j\).
• Collocation (bem): test function = Dirac delta at the segment midpoint. The reaction matrix entry \(Z_{ij}\) is the point-evaluated potential of segment \(i\) at its midpoint caused by a unit current on segment \(j\).

In the grounding literature, Colominas, Navarrina & Casteleiro (2007, 2012) document the collocation BEM as the historically dominant variant for layered-soil grounding analysis. It has roughly half the cost of Galerkin per matrix entry (one integration instead of a double integration) and gives essentially identical accuracy on smooth electrodes. The price is slightly higher sensitivity to the segment-length / wire-radius ratio at the wire end-points.

bem was added to the engine family to provide a methodological alternative to the Galerkin scheme. That intent is only partly realised: on the diagonal both engines use the same analytical line self-potential and off the diagonal both reduce to a point-source evaluation at the segment midpoint, so on the soils bem accepts (\(n \le 2\)) the two assembled matrices coincide to floating-point noise. The distinction between Galerkin and collocation is real in theory but has no numerical footprint in this implementation — see the status box above.

Governing equation: boundary integral

The same boundary integral equation as in mom:

\[ \sum_{j=1}^{N} Z_{ij}\, I_j \;=\; \varphi_c \qquad \forall\, i \in c, \qquad \sum_{j \in c} I_j \;=\; I_{c,\text{in}}, \]

with the reaction matrix entries

\[ Z_{ij} \;=\; \frac{1}{4\pi} \int_{\Sigma_j} G(\mathbf{r}_i, \mathbf{r}'_j)\, dS_j. \]

The difference: \(\mathbf{r}_i\) is now the centre of segment \(i\) rather than the average over its length. The \(j\)-side integration is unchanged.

For thin-wire grounding electrodes (radius \(a \ll L_i\)) the \(\Sigma_j\) surface integral collapses to a line integral with appropriate kernel; for the off-diagonal entries the line integral itself further reduces to a point-source evaluation at the segment midpoint (the segment is short compared to its distance to the field point). The diagonal carries the analytical line self-potential, the same as in mom and image.

Numerical strategy

Kernel choice

bem uses the same kernel infrastructure as the rest of the engine family. For each soil class:

• HomogeneousSoil → homogeneous self-kernel (\(1/r + 1/r_{\text{air}}\) point-source off-diagonal, line self-potential on the diagonal). Bit-exact match to image and mom at the Galerkin level for \(n = 1\).
• TwoLayerSoil → the closed-form Tagg / Sunde self-kernel (_two_layer_self_kernel_factory, with allow_cross_layer=True so interface-crossing geometries take the rigorous ADR-0007 path). At \(n = 2\) the constant \(\Gamma_1 \equiv K_1\) is a single complex image at \(\beta = 0\), so the geometric series is the exact complex-image representation — there is nothing for a fit to add.
• MultiLayerSoil (\(n \ge 3\)) → NotImplementedError. The historic contribution for this regime came from an incomplete Green's function (no \(2 h_1\) image families, no surface-interface multiple-reflection denominator; see solver/_layered.py). Use mom_sommerfeld or fem.

No complex-image fit is computed. Up to 0.14.1 solve_bem called fit_complex_images on every solve and published its (failed) diagnostics as cim_n_images / cim_rms, although the only reachable branches never looked at them (review pass 9, F34). The metadata now records cim_fit_used = False, cim_n_images = 0, cim_rms = None and reduces_to = "mom".

The end result: bem and mom differ only in the test function, and on the soils bem accepts even that difference cancels — the matrices are identical.

Reaction matrix assembly

For the homogeneous and 2-layer cases the assembly is one call to the existing self-kernel factory with the identity matrix as the "currents" argument — the resulting matrix is exactly \(Z\). The diagonal carries the line self-potential; the off-diagonals carry the point-source approximation.

The \(n \ge 3\) branch used to add a complex-image contribution

\[ Z^{(\text{layered})}_{ij} \;=\; Z^{(\text{hom})}_{ij} + \frac{\rho_1}{4\pi}\, \sum_{k=1}^{P} \frac{a_k}{\sqrt{s_{ij}^2 + (z_i + z_j + 2\beta_k)^2}}, \]

with \(s_{ij}\) the radial distance between the segment midpoints and \(a_k, \beta_k\) the matrix-pencil fit coefficients. That branch has been unreachable since 0.11.0 (the \(n \ge 3\) rejection sits in front of it) and was deleted in 0.15.0: the sum represents only a single image family of the layered Green's function, so keeping it alive invited a silent wrong-physics path. The matrix builder now raises for \(n \ge 3\) instead. The formula stays documented here because it is the shape a complete \(n \ge 3\) kernel would take once the missing families are added.

Linear-system solve

The cluster augmenting rows and the Galerkin solve are reused from mom — the only difference between mom and bem in the assembled system is the matrix entries themselves.

Validity envelope

Property	Range / value
Soil model	HomogeneousSoil, TwoLayerSoil (a MultiLayerSoil is accepted only while it reduces to \(n \le 2\))
Frequency	quasi-static, \(f < 1\,\text{kHz}\)
Number of layers	\(n \le 2\); \(n \ge 3\) raises NotImplementedError
Electrode placement	free — \(n = 2\) interface crossings dispatch to the ADR-0007 cross-layer kernel
Wire radius / segment ratio	thin-wire, \(a \ll L_i\)
Mesh size \(N\)	\(\le 1000\) at acceptable runtime
Number of complex images \(P\)	not applicable — no fit runs

Convergence and cost

• Per-segment accuracy. Equal to mom on the soils this backend accepts: the assembled matrices agree to 4e-16, so the two engines converge along the segment-length axis in lockstep. (The textbook difference between collocation and Galerkin — slightly faster convergence, slightly higher sensitivity to the wire-radius / segment-length ratio at wire ends — would only appear with a genuinely averaged Galerkin kernel.)
• Computational cost. \(O(N^2)\) matrix build, \(O((N + K)^3)\) solve.
• Reduction. At \(K_1 = 0\) the engine collapses bit-exactly to the homogeneous bem solution, which itself agrees with image to within the segment-discretisation envelope.

Cross-validation notes

Counterpart	Expected agreement	What is checked
mom (\(n \le 2\))	bit-exact (4e-16)	same computation — identical reaction matrix and constraint solve; no independent information
image (\(n = 1\))	\(\le 5\,\%\)	uniform-current vs. collocation weighting
image_2layer (\(n = 2\))	\(\le 5\,\%\)	uniform-current vs. collocation weighting on the same Tagg / Sunde kernel
mom_sommerfeld (\(n \le 2\))	\(\le 5\,\%\)	direct quadrature of the full layered Green's function — the genuinely independent kernel
fem (\(n \le 2\))	\(\le 10\,\%\)	volume PDE, independent problem form
Sunde / Dwight closed forms	\(\le 5\,\%\)	tighter than the image backend

bem is a flavour of the integral-equation family, not an independent line of defence: the historic "mom / bem / cim triangle" collapses to a single point, because bem reproduces mom bit-for-bit and cim reproduces image_2layer bit-for-bit. The independent checks are mom_sommerfeld (quadrature of the full layered kernel) and fem (volume PDE); see the ADR-0002 amendment.

References

• Colominas, I., Navarrina, F. & Casteleiro, M. (2007). Numerical simulation of transferred potentials in earthing grids considering layered soil models. IEEE PWRD 22(3). Layered BEM for grounding systems.
• Colominas, I., París, J., Navarrina, F. & Casteleiro, M. (2012). Improvement of computer methods for grounding analysis in layered soils by using high-efficient convergence acceleration techniques. Adv. Eng. Soft. 44. Aitken / Pade acceleration of the BEM kernel; cross-checks against measurement.
• Brebbia, C. A. & Dominguez, J. (1992). Boundary Elements: An Introductory Course, McGraw-Hill. The BEM textbook.
• Harrington, R. F. (1968). Field Computation by Moment Methods, Macmillan. Cross-reference for the Galerkin alternative.

Example

import groundfield as gf

soil = gf.TwoLayerSoil(rho_1=100.0, rho_2=500.0, h_1=2.0)
world = gf.create_world(soil=soil)
gf.create_electrode(world, "rod", name="g1",
                    position=(0.0, 0.0, 0.0), length=1.5)
gf.create_source(world, attached_to="g1", magnitude=1.0)

engine = gf.create_engine(backend="bem",
                          segment_length=0.1,
                          frequencies=[50.0])
result = world.solve(engine)
print(result.cluster_impedance("g1")[0])
print(result.metadata["cim_fit_used"])   # False -- no fit runs
print(result.metadata["reduces_to"])     # 'mom'

For \(n \ge 3\) the call raises; use mom_sommerfeld (full layered Green's function) or fem (volume PDE) instead.

API reference

bem

Boundary-Element-Method backend (bem).

Mathematical / physical model

Following Colominas, Navarrina & Casteleiro (2007, 2012) we treat the grounding system as a boundary problem: the unknown is the leakage-current density along the wire surfaces. Discretising every electrode into \(N\) line segments turns the boundary integral equation into a dense linear system $$ \sum_{j=1}^{N}\, Z_{ij}\, I_j \;=\; \varphi_c \qquad \forall\,i \in c, \qquad \sum_{j \in c} I_j \;=\; I_{c,\text{in}}, $$ where \(c\) ranges over the galvanic clusters, \(\varphi_c\) is the (unknown) shared cluster potential, and the reaction matrix entries $$ Z_{ij} \;=\; \frac{1}{4\pi} \int_{\Sigma_i} G(\mathbf{r}_i, \mathbf{r}'_j)\, dS_j $$ are obtained by collocation of \(G\) at the centre of segment \(i\) (single test point per segment), instead of the Galerkin average used in :mod:groundfield.solver.mom. Collocation is the historically dominant flavour of BEM in the grounding literature because it preserves the same accuracy on smooth electrodes while being roughly half the cost of the Galerkin scheme.

The Green's function \(G\) is the closed-form self-kernel of the image-charge family: \(G(r) = 1/r + 1/r_{\text{air-img}}\) for homogeneous soil, and the exact Tagg/Sunde series of :mod:groundfield.solver.image_2layer for a two-layer soil.

Differences from mom — and why they are nil in practice

(Corrected in the 2026-07-08 audit, WP-E/P4; the last remnant of the complex-image story removed in 0.15.0, review pass 9 F34.) For n <= 2 — the only regime this backend accepts — the two engines assemble the identical reaction matrix from the same kernels and solve it with the same constraint solver; the measured relative difference is 4e-16, i.e. bem reproduces mom bit-for-bit up to floating-point associativity. Their mutual agreement therefore validates the shared discretisation, not the physics, and a cross-validation table must not count them as two engines. Genuine methodological independence in the layered cross-checks comes from mom_sommerfeld (direct Sommerfeld quadrature) and fem (volume PDE). bem is kept as the collocation-flavoured entry point of the family. n >= 3 soils are rejected — the historic complex-image kernel this backend shared with cim was structurally incomplete — and consequently no complex-image fit is computed any more: up to 0.14.1 solve_bem called fit_complex_images on every solve and published its (failed) diagnostics as cim_n_images / cim_rms although no reachable path consumed them.

Validity

• Quasi-static, \(f < 1\,\mathrm{kHz}\).
• HomogeneousSoil and TwoLayerSoil only (MultiLayerSoil is accepted while it reduces to \(n \le 2\)).
• Wire radius small compared to segment length; thin-wire approximation in the line self-correction (same as the other segment-based engines).

References

• Colominas, I., Navarrina, F., & Casteleiro, M. (2007). Numerical simulation of transferred potentials in earthing grids considering layered soil models. IEEE PWRD 22(3).
• Colominas, I., París, J., Navarrina, F., & Casteleiro, M. (2012). Improvement of computer methods for grounding analysis in layered soils by using high-efficient convergence acceleration techniques. Adv. Eng. Soft. 44.

solve_bem

solve_bem(world: 'World', engine: 'Engine') -> FieldResult

Boundary-Element-Method solver (collocation, closed-form kernels).

What actually runs

• \(n = 1\) → homogeneous image-charge kernel.
• \(n = 2\) → exact Tagg/Sunde series kernel.
• \(n \ge 3\) → :class:NotImplementedError (audit 2026-07-08, WP-E). Use mom_sommerfeld or fem.

In both reachable regimes the assembled reaction matrix and the constraint solve are the same as in mom (relative difference 4e-16), so this result is not an independent check of mom — see the module docstring and the ADR-0002 amendment. No complex-image fit is computed (review pass 9, F34).

Parameters:

Name	Type	Description	Default
world	'World'	World to evaluate.	required
engine	'Engine'	Engine configuration; engine.segment_length controls the discretisation, engine.image_max_terms / engine.image_series_tol the Tagg/Sunde truncation at \(n = 2\).	required

Returns:

Type	Description
FieldResult	metadata['cim_fit_used'] = False and metadata['reduces_to'] = 'mom' record both facts above.

Source code in src/groundfield/solver/bem.py

def solve_bem(world: "World", engine: "Engine") -> FieldResult:
    """Boundary-Element-Method solver (collocation, closed-form kernels).

    What actually runs
    ------------------
    - $n = 1$ → homogeneous image-charge kernel.
    - $n = 2$ → exact Tagg/Sunde series kernel.
    - $n \\ge 3$ → :class:`NotImplementedError` (audit 2026-07-08,
      WP-E). Use ``mom_sommerfeld`` or ``fem``.

    In both reachable regimes the assembled reaction matrix and the
    constraint solve are the same as in ``mom`` (relative difference
    4e-16), so this result is **not** an independent check of ``mom``
    — see the module docstring and the ADR-0002 amendment. No
    complex-image fit is computed (review pass 9, F34).

    Parameters
    ----------
    world
        World to evaluate.
    engine
        Engine configuration; ``engine.segment_length`` controls the
        discretisation, ``engine.image_max_terms`` /
        ``engine.image_series_tol`` the Tagg/Sunde truncation at
        $n = 2$.

    Returns
    -------
    FieldResult
        ``metadata['cim_fit_used'] = False`` and
        ``metadata['reduces_to'] = 'mom'`` record both facts above.
    """
    if not isinstance(world.soil, (HomogeneousSoil, TwoLayerSoil, MultiLayerSoil)):
        raise TypeError(
            "Backend 'bem' supports HomogeneousSoil, TwoLayerSoil, "
            f"and MultiLayerSoil. Got: {type(world.soil).__name__}."
        )
    if not world.electrodes:
        raise ValueError("World contains no electrodes.")
    _reject_concrete_shells(world, "bem")
    _warn_ignored_sources(world, "bem")

    stack = as_layer_stack(world.soil)
    if stack.n_layers >= 3:
        # Audit 2026-07-08, WP-E: bem shares the historic incomplete
        # n>=3 CIM kernel — reject loudly until a complete kernel
        # exists.
        raise NotImplementedError(
            f"bem: n_layers = {stack.n_layers} >= 3 is not supported — "
            "the shared complex-image kernel was structurally "
            "incomplete (audit 2026-07-08). Use "
            "backend='mom_sommerfeld' (full layered Green's function) "
            "or 'fem' for n >= 3 soils."
        )
    ds = engine.segment_length

    _log.info(
        "bem: n_layers=%d, segment_length=%.3f, closed-form kernel "
        "(no complex-image fit); identical reaction matrix to 'mom'",
        stack.n_layers, ds,
    )

    # 1) Discretisation.
    all_segments: list[_Segment] = []
    elec_to_segidx: dict[str, list[int]] = {}
    interfaces = (
        (float(stack.h[0]),) if stack.n_layers >= 2 else None
    )
    for e in world.electrodes:
        segs = _discretize_electrode(e, ds, layer_interfaces=interfaces)
        elec_to_segidx[e.name] = list(
            range(len(all_segments), len(all_segments) + len(segs))
        )
        all_segments.extend(segs)

    # 2) Per-electrode input currents.
    elec_input_current: dict[str, complex] = {
        e.name: 0j for e in world.electrodes
    }
    for src in world.sources:
        if src.kind != "current":
            continue
        i_complex = src.magnitude * np.exp(1j * np.deg2rad(src.phase_deg))
        if src.attached_to in elec_input_current:
            elec_input_current[src.attached_to] += i_complex

    cluster_id = _build_clusters(world.electrodes, world.conductors)
    finite_branches = _build_finite_branches(world.conductors, cluster_id)

    # 2b) Distributed-conductor topology (ADR-0003) + ADR-0004
    #     inductive coupling.
    cond_segs, distributed_branches_objs, interior_nodes = _build_distributed_topology(
        world.conductors, cluster_id
    )
    for s in cond_segs:
        pn = s.electrode_name
        elec_to_segidx[pn] = [len(all_segments)]
        all_segments.append(s)
        cluster_id[pn] = pn
    for n_ in interior_nodes:
        if n_ not in cluster_id:
            cluster_id[n_] = n_
            elec_to_segidx[n_] = []
    n_lumped_branches = len(finite_branches)
    distributed_branch_tuples = [
        (db.node_a, db.node_b, db.R) for db in distributed_branches_objs
    ]
    finite_branches = list(finite_branches) + distributed_branch_tuples
    earth_inductive_model = getattr(
        engine, "earth_inductive_model", "perfect_mirror"
    )
    sigma_earth_for_carson: float | None = None
    layered_earth_for_sommerfeld: object = None
    if earth_inductive_model == "carson_series":
        from groundfield.coupling import resolve_earth_conductivity

        sigma_earth_for_carson = resolve_earth_conductivity(world.soil)
    elif earth_inductive_model == "sommerfeld":
        from groundfield.coupling import resolve_earth_layers

        layered_earth_for_sommerfeld = resolve_earth_layers(world.soil)
    inductance_matrix_full, has_inductance, carson_builder = _assemble_inductance_matrix(
        distributed_branches_objs,
        n_lumped_branches=n_lumped_branches,
        n_total_branches=len(finite_branches),
        earth_model=earth_inductive_model,
        sigma_earth=sigma_earth_for_carson,
        layered_earth=layered_earth_for_sommerfeld,
    )

    n_segments = len(all_segments)
    seg_points = np.array([s.midpoint for s in all_segments])
    seg_lengths = np.array([s.length for s in all_segments])
    wire_radii = np.array([s.wire_radius for s in all_segments])

    # 3) Reaction matrix via collocation (series knobs from engine).
    #    The n>=3 cross-layer guard that used to sit here was
    #    unreachable behind the n>=3 rejection above and was removed in
    #    0.15.0; n=2 cross-layer geometries are handled rigorously by
    #    _two_layer_self_kernel_factory(allow_cross_layer=True)
    #    (ADR-0007).
    Z = _build_Z_collocation(
        seg_points, seg_lengths, wire_radii, stack,
        max_terms=engine.image_max_terms, tol=engine.image_series_tol,
    )

    # 4) Frequency loop (Galerkin solve + Z · I_seg for phi).
    n_freq = len(engine.frequencies)
    omegas = [2.0 * np.pi * float(f) for f in engine.frequencies]
    real_electrode_names = {e.name for e in world.electrodes}

    def _solve_at(omega: float) -> tuple[np.ndarray, np.ndarray]:
        carson_dz = (
            carson_builder(omega) if (has_inductance and carson_builder is not None)
            else None
        )
        sc, _ = _galerkin_solve(
            Z=Z,
            elec_input_current=elec_input_current,
            cluster_id=cluster_id,
            elec_to_segidx=elec_to_segidx,
            n_segments=n_segments,
            finite_branches=finite_branches,
            omega=omega if has_inductance else 0.0,
            inductance_matrix=inductance_matrix_full if has_inductance else None,
            carson_correction=carson_dz,
        )
        ph = np.zeros(n_segments, dtype=complex)
        if sc.any():
            ph = Z @ sc.real + 1j * (Z @ sc.imag)
        return sc, ph

    sc_per_freq: list[np.ndarray] = []
    phi_per_freq: list[np.ndarray] = []
    if has_inductance:
        for omega in omegas:
            sc, ph = _solve_at(omega)
            sc_per_freq.append(sc)
            phi_per_freq.append(ph)
    else:
        sc, ph = _solve_at(0.0)
        sc_per_freq = [sc] * n_freq
        phi_per_freq = [ph] * n_freq

    electrode_potentials: dict[str, list[complex]] = {}
    electrode_currents: dict[str, list[complex]] = {}
    conductor_currents: dict[str, list[complex]] = {}
    conductor_potentials: dict[str, list[complex]] = {}
    for ename, idxs in elec_to_segidx.items():
        if not idxs:
            continue
        u_list = [
            complex(np.mean(phi_per_freq[k][idxs])) for k in range(n_freq)
        ]
        i_list = [
            complex(sc_per_freq[k][idxs].sum()) for k in range(n_freq)
        ]
        if ename in real_electrode_names:
            electrode_potentials[ename] = u_list
            electrode_currents[ename] = i_list
        else:
            conductor_potentials[ename] = u_list
            conductor_currents[ename] = i_list

    point_sources = [
        PointSource(
            position=tuple(seg_points[i].tolist()),
            current=[complex(sc_per_freq[k][i]) for k in range(n_freq)],
            electrode_name=all_segments[i].electrode_name,
            length=float(seg_lengths[i]),
        )
        for i in range(n_segments)
    ]
    cluster_members: dict[str, list[str]] = {}
    for ename in real_electrode_names:
        cluster_members[ename] = sorted(
            n for n in cluster_id
            if cluster_id[n] == cluster_id[ename] and n in real_electrode_names
        )

    metadata = {
        "world_name": world.name,
        "n_segments": n_segments,
        "segment_length": ds,
        "n_layers": int(stack.n_layers),
        "rhos": stack.rhos.tolist(),
        "h": stack.h.tolist(),
        # Honesty flags (review pass 9, F34): the reachable n <= 2
        # paths use exact closed-form kernels, so there is no
        # complex-image fit to report — and the reaction matrix is the
        # one `mom` assembles, so this is not an independent check.
        "cim_fit_used": False,
        "cim_n_images": 0,
        "cim_rms": None,
        "reduces_to": "mom",
        "solver": "collocation",
        "stub": False,
        "earth_inductive_model": earth_inductive_model,
    }
    if has_inductance:
        from groundfield.coupling.carson import skin_depth

        sigma_ref = (
            sigma_earth_for_carson
            if sigma_earth_for_carson is not None
            else 1.0 / float(stack.rhos[0])
        )
        metadata["penetration_depth"] = {
            float(f): skin_depth(2.0 * np.pi * f, sigma_ref)
            for f in engine.frequencies
        }
    if conductor_currents:
        metadata["conductor_node_currents"] = conductor_currents
        metadata["conductor_node_potentials"] = conductor_potentials

    return FieldResult(
        backend="bem",
        frequencies=list(engine.frequencies),
        electrode_potentials=electrode_potentials,
        electrode_currents=electrode_currents,
        point_sources=point_sources,
        soil_resistivity=float(stack.rhos[0]),
        soil=world.soil,
        clusters=cluster_members,
        metadata=metadata,
    )

Related material

• ADR-0002 — engine selection heuristic, and the 2026-07-09 amendment that revoked bem's role as an independent cross-validation engine.