bem — Boundary-Element Collocation

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. When mom and bem agree on a given problem, the answer is robust against the choice of test function; when they disagree, the disagreement is reproducible and quantifiable.

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). The matrix-pencil-fit-based CIM kernel is intentionally not used here, even though bem takes a fit_complex_images instance: at \(n = 2\) the \(\Gamma_1 \equiv K_1\) constant makes the matrix-pencil fit ill-conditioned (a single pole at \(\beta = 0\)). Falling back on the geometric series gives a bit-exact match with image_2layer.
• MultiLayerSoil (\(n \ge 3\)) → homogeneous self-kernel for the direct + air-mirror part plus the closed-form complex-image contribution from cim. This is the same closed-form layered Green's function the cim engine uses, but evaluated with collocation rather than Galerkin averaging.

The end result: bem and mom differ only in the test function; the kernel is identical for \(n = 1, 2\), and shares the same CIM approximation for \(n \ge 3\).

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.

For \(n \ge 3\) the homogeneous part is built first, then the complex-image contribution is added explicitly:

\[ 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. The imaginary part of the sum cancels by symmetry, so we take the real part to suppress numerical residue.

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, MultiLayerSoil
Frequency	quasi-static, \(f < 1\,\text{kHz}\)
Electrode placement	every segment in the upper layer
Wire radius / segment ratio	thin-wire, \(a \ll L_i\)
Mesh size \(N\)	\(\le 1000\) at acceptable runtime
Number of complex images \(P\)	inherited from cim (default 8)

Convergence and cost

• Per-segment accuracy. Comparable to mom for smooth electrodes (rods, rings, meshes); collocation converges somewhat faster on the segment-length axis but is slightly more sensitive to the wire-radius / segment-length ratio at the electrode ends. For the typical geometries the difference is negligible.
• Computational cost. \(O(N^2)\) matrix build, \(O((N + K)^3)\) solve. For the layered case the matrix build is dominated by the one-shot CIM fit (\(O(N_s P^2)\), \(N_s = 64\), \(P = 8\), so negligible).
• 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
image (\(n = 1\))	\(\le 5\,\%\)	uniform-current vs. collocation
image_2layer (\(n = 2\))	bit-exact	bem reuses the Tagg / Sunde kernel for \(n = 2\)
mom (any \(n\))	\(\le 5\,\%\)	Galerkin vs. collocation on the same kernel
cim (\(n \ge 3\))	\(\le 5\,\%\)	shares the CIM kernel; only the test function differs
mom_sommerfeld (any \(n\))	\(\le 5\,\%\)	quadrature reference
Sunde / Dwight closed forms	\(\le 5\,\%\)	tighter than the image backend

bem is the alternative-weighting cross-check in the matrix. It is paired with mom and cim in the test suite; together they form a triangle that detects bugs in the kernel (mom vs. cim disagreement), the test function (mom vs. bem disagreement), or the layered approximation (cim vs. mom_sommerfeld disagreement).

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])

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 layered Sommerfeld kernel itself; we evaluate it through the closed-form complex-image fit provided by :mod:groundfield.solver.cim, so the BEM and the CIM engines share the same physics for layered soils. For homogeneous soils \(G(r) = 1/r + 1/r_{\text{air-img}}\).

Differences from mom

(Corrected in the 2026-07-08 audit, WP-E/P4.) For n <= 2 the two engines assemble the identical reaction matrix from the same kernels and solve it with the same constraint solver — measured relative difference 4e-16. Their mutual agreement therefore validates the shared discretisation, not the physics; 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 CIM kernel was structurally incomplete).

Validity

• Quasi-static, \(f < 1\,\mathrm{kHz}\).
• 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",
    *,
    n_images: int = 8,
    n_samples: int = 64
) -> FieldResult

Boundary-Element-Method solver (collocation, layered CIM kernel).

Parameters:

Name	Type	Description	Default
world	'World'	World to evaluate.	required
engine	'Engine'	Engine configuration.	required
n_images	int	Forwarded to :func:groundfield.solver.cim.fit_complex_images when the soil is layered. Ignored otherwise.	8
n_samples	int	Forwarded to :func:groundfield.solver.cim.fit_complex_images when the soil is layered. Ignored otherwise.	8

Returns:

Type	Description
FieldResult

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

def solve_bem(
    world: "World",
    engine: "Engine",
    *,
    n_images: int = 8,
    n_samples: int = 64,
) -> FieldResult:
    """Boundary-Element-Method solver (collocation, layered CIM kernel).

    Parameters
    ----------
    world
        World to evaluate.
    engine
        Engine configuration.
    n_images, n_samples
        Forwarded to :func:`groundfield.solver.cim.fit_complex_images`
        when the soil is layered. Ignored otherwise.

    Returns
    -------
    FieldResult
    """
    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."
        )
    fit = fit_complex_images(stack, n_images=n_images, n_samples=n_samples)
    ds = engine.segment_length

    _log.info(
        "bem: n_layers=%d, n_images=%d, segment_length=%.3f",
        stack.n_layers, fit.a.size, 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])

    if stack.n_layers >= 3:
        z_max = seg_points[:, 2].max()
        h_1 = float(stack.h[0])
        if z_max >= h_1:
            # n>=3 cross-layer is not implemented — the complex-image
            # kernel would silently apply the upper-layer expansion
            # below the first interface. Hard error instead of a
            # warning followed by wrong physics. (n=2 cross-layer is
            # handled via the _two_layer_self_kernel_factory
            # dispatcher with allow_cross_layer=True.)
            raise ValueError(
                f"bem: cross-layer geometry on n_layers="
                f"{stack.n_layers} is not supported (z_max = "
                f"{z_max:.3f} m >= h_1 = {h_1:.3f} m). Use "
                "backend='image_2layer' for n=2 cross-layer worlds; "
                "for n>=3 use 'mom_sommerfeld' or thicken the upper "
                "layer."
            )

    # 3) Reaction matrix via collocation (series knobs from engine).
    Z = _build_Z_collocation(
        seg_points, seg_lengths, wire_radii, stack, fit,
        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(),
        "cim_n_images": int(fit.a.size),
        "cim_rms": float(fit.rms),
        "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; bem is the alternative-weighting cross-check in the layered family.