Local Domain Boundary Element Method

• LD-BEM is a boundary-integral technique that partitions a global problem into local subdomains, enforcing governing equations via boundary operators.
• The method employs local elimination of interior fields, yielding sparse, efficient global systems ideal for time-fractional and convection–diffusion–reaction problems.
• LD-BEM extends classical BEM by handling nonlinearities, memory terms, and interface coupling, making it versatile for applications like electromagnetic scattering and metasurface analysis.

Local Domain Boundary Element Method (LD-BEM) denotes a class of boundary-integral formulations in which a global problem is decomposed into local subdomains or elements, the governing equations are enforced by boundary operators or local-domain integral equations on each part, and the global solution is recovered by coupling only interface or skeletal unknowns. In this sense, LD-BEM is not a single algebraic template but a boundary-centered computational paradigm that appears in nonlinear time-fractional diffusion–reaction, convection–diffusion–reaction boundary value problems, and electromagnetic scattering in interconnected homogeneous regions with metasurface interfaces (Gortsas, 22 Aug 2025, Hofreither et al., 2015, Stewart et al., 2018).

1. Conceptual basis and departure from conventional BEM

The classical Boundary Element Method is especially attractive when a fundamental solution is available and the governing problem can be reduced to boundary-only integral equations. For diffusion–reaction equations this is the case for the Laplacian in the absence of domain sources, but the situation changes decisively for nonlinear and time-fractional models: there is no fundamental solution for the nonlinear operator, so the operator is split into a linear diffusion part and a nonlinear source term, while the fractional time derivative introduces a nonlocal-in-time contribution that, under a static fundamental-solution approach, appears as a domain source at each time step. The result is the familiar dense BEM algebra with quadratic storage and cubic direct-solver cost, together with repeated evaluation of volume integrals and nonlinear iterations (Gortsas, 22 Aug 2025).

LD-BEM modifies this architecture by moving from a single global boundary representation to a collection of local boundary representations. In the two-dimensional time-fractional Fisher–KPP setting, the domain is partitioned into conformal, non-overlapping subregions, the boundary integral equation is written on each subregion separately, and the domain integrals are converted into local-domain integrals attached to each local cell. In the three-dimensional convection–diffusion–reaction setting, the same principle appears in a skeletal formulation: if $L u = 0$ in an element $T$, then

$$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$

so the volumetric form is replaced by a boundary pairing with a local Dirichlet-to-Neumann, or Steklov–Poincaré, operator $S_T$ (Hofreither et al., 2015).

A common misconception is that boundary-element formulations remain boundary-only even after nonlinearities, memory terms, or complicated interface physics are introduced. The local-domain literature shows the opposite. When domain contributions are unavoidable, LD-BEM reorganizes them so that they remain local and are eliminated or condensed at the subregion level rather than propagated into a globally dense coupling.

2. Local partitioning, interface coupling, and algebraic structure

Across the formulations documented in the cited works, LD-BEM is characterized by local elimination of interior fields or values and global coupling through interfaces, skeleton traces, or shared current unknowns (Gortsas, 22 Aug 2025, Hofreither et al., 2015, Stewart et al., 2018).

Setting	Local object	Global coupling
Time-fractional Fisher–KPP	Quadrilateral subregion with boundary and local-domain matrices $[H]$, $[G]$, $[C]$	Sparse, non-symmetric assembled system on subregion boundaries
Convection–diffusion–reaction	Polyhedral element with local Steklov–Poincaré operator $S_T$	Sparse skeleton system assembled from $$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$
Electromagnetic metasurfaces	Homogeneous region represented by boundary current operators $\mathcal{L}$, $T$0	Dense interface-current system coupled by GSTCs or conventional interface conditions

In the time-fractional Fisher–KPP formulation, interface conditions between neighboring subregions are explicit: $T$1 Unknowns are restricted to boundary degrees of freedom on subregion interfaces and outer boundaries, while interior values are eliminated at the subregion level. In the reported implementation, each quadrilateral cell has four interior nodes and eight boundary nodes because linear discontinuous boundary elements are used on the four edges. This produces local $T$2 boundary-to-boundary blocks, $T$3 boundary-to-interior blocks, and $T$4 and $T$5 domain blocks before local elimination.

In the convection-adapted BEM-based FEM, the coupling variable is not a flux/value pair on explicit cell boundaries but the trace on the mesh skeleton $T$6. The discrete skeletal space $T$7 is spanned by nodal basis functions that are PDE-harmonic inside each element and whose traces on edges and faces solve lower-dimensional projected PDEs. In the metasurface scattering formulation, each local homogeneous region is represented only by surface currents on its boundary, and coupling occurs on shared interfaces through Generalized Sheet Transition Conditions (GSTCs) or through conventional PEC, PMC, or dielectric interface relations.

This suggests that LD-BEM is best understood as a unifying principle—local boundary reduction plus interface coupling—rather than as a single fixed discretization recipe.

3. Time-fractional nonlinear diffusion–reaction formulations

A detailed contemporary realization of LD-BEM is the numerical treatment of the two-dimensional time-fractional nonlinear Fisher–KPP equation,

$T$8

with

$T$9

supplemented by Dirichlet, Neumann, and, where appropriate, Robin boundary conditions. The formulation accommodates the Caputo derivative and the fractal-fractional derivative in the Riemann–Liouville sense for $$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$0. The spatial representation uses the static Laplace fundamental solution

$$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$1

and singular and near-singular boundary integrals are evaluated with the semi-analytical procedures of Guiggiani and Gigante. The Caputo derivative is discretized by a first-order accurate L1-type scheme on a uniform grid,

$$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$2

with a convolution history term

$$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$3

The nonlinear logistic term is treated by a lagging fixed-point linearization,

$$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$4

so only the iteration-dependent blocks associated with the nonlinear reaction must be updated, while the geometry-dependent $$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$5, $$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$6, and $$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$7 matrices can be reused (Gortsas, 22 Aug 2025).

The algorithmic structure is local and sparse. At each time step, the fractional history term is updated, local matrices are assembled on each quadrilateral subregion, the reaction term is linearized by lagging, interior unknowns are eliminated by inverting $$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$8 with an iterative Sherman–Morrison-type expansion, and a sparse, non-symmetric global system

$$a_T(u,v)=\langle \gamma_T^1 u,\gamma_T^0 v\rangle_{\partial T} =\langle S_T\,u_{\partial T},v_{\partial T}\rangle_{\partial T},$$9

is assembled and solved for boundary unknowns. Interior values are then recovered by back-substitution. The history term grows linearly with the number of time steps because it is implemented explicitly and without compression.

The reported numerical evidence covers six two-dimensional problems. On the rectangular manufactured-solution test with $S_T$0, $S_T$1 quadrilateral subregions, $S_T$2, and $S_T$3–$S_T$4, LD-BEM with constant and linear local approximations matches the analytical solution well; for the linear case at $S_T$5, the reported errors are $S_T$6 and $S_T$7. In the solid-cylinder radial diffusion problem, LD-BEM and FPM agree very well with the analytical solution at $S_T$8 for $S_T$9; one reported LD-BEM value is $[H]$0 at $[H]$1. Additional tests on a hollow cylinder, a large rectangle, an intersection-of-circles domain, and a star-shaped domain show agreement with analytical solutions or with an independent meshless FPM, together with stable nonlinear iterations for the tested parameter ranges. Across these cases, decreasing $[H]$2 slows diffusion and front propagation at fixed physical time, consistent with fractional-memory induced subdiffusion.

4. Skeletal and Trefftz interpretations in convection–diffusion–reaction

In three-dimensional convection–diffusion–reaction problems, LD-BEM appears in a distinctly skeletal form. The governing operator is

$[H]$3

with $[H]$4 symmetric and uniformly positive definite, $[H]$5, $[H]$6, and $[H]$7 in $[H]$8, $[H]$9, and piecewise-constant coefficients on a polyhedral mesh. The method is described as a non-standard finite element method with PDE-harmonic shape functions on polyhedral elements, constructed by means of local boundary element techniques; it can therefore be considered a local Trefftz method with element-wise PDE-harmonic shape functions. The local Steklov–Poincaré operator is represented as

$[G]$0

and the element stiffness contribution is

$[G]$1

The traces of the basis functions are not piecewise-linear polynomials. Along edges they solve a projected one-dimensional operator $[G]$2, and on faces they solve a projected two-dimensional convection–diffusion–reaction problem $[G]$3. Because faces may be convection-dominated, the construction uses SUPG stabilization on auxiliary face triangulations, with the local stabilization parameter chosen from a Peclet-number criterion: $[G]$4 when $[G]$5, and $[G]$6 otherwise. The approximate element boundary stiffness matrix is

$[G]$7

and the global system is assembled as

$[G]$8

The resulting linear system is sparse and non-symmetric and is solved by GMRES (Hofreither et al., 2015).

The method’s principal significance lies in convection adaptation. Earlier BEM-based FEM variants used piecewise-linear face traces, whereas the convection-adapted scheme constructs traces from projected edge and face problems aligned with the local convection field and reinforced by SUPG stabilization and Shishkin-like face-mesh adaptation. The rationale parallels residual-free bubbles: with exact DtN evaluation, the method coincides with RFB FEM in diffusion-only settings, and the Trefftz character suppresses interior residuals because basis functions satisfy $[G]$9 elementwise.

The numerical evidence is specifically framed in terms of stability under increasing element Peclet number. On a tetrahedral mesh of the unit cube with $[C]$0 elements, $[C]$1 faces, $[C]$2 edges, and $[C]$3 nodes, standard FEM is reported stable up to $[C]$4, the previous BEM-based FEM up to $[C]$5, and the new convection-adapted method up to $[C]$6 for $[C]$7 and $[C]$8 for $[C]$9. On the same problem, GMRES iteration counts remain lower and more stable than in the previous scheme; at $S_T$0, $S_T$1 iterations are reported for $S_T$2 versus $S_T$3 for the previous method. On a polyhedral prismatic mesh with a rotating convection field, stability deteriorates earlier because of complex face geometries, but geometric row scaling yields moderate and controlled GMRES counts across the Peclet range.

5. Interconnected-region electromagnetic formulations and metasurface interfaces

In frequency-domain electromagnetics, LD-BEM is realized as a multi-region boundary-current formulation for homogeneous subdomains coupled only through shared interfaces. For a region with wavenumber $S_T$4, the scattered fields are represented by equivalent electric and magnetic surface currents $S_T$5 and $S_T$6: $S_T$7 where $S_T$8 and $S_T$9 are boundary integral operators built from the Helmholtz Green’s function. The formulation is given for both 3D and 2D kernels, although the reported demonstrations are in 2D. Metasurfaces are modeled as zero-thickness boundaries with effective sheet susceptibilities and enforced through GSTCs. In scalar mono-isotropic form,

$$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$0

with diagonal matrices $$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$1 and $$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$2 containing the elementwise sheet parameters. In two dimensions, interfaces are discretized into line segments with pulse basis functions and collocation at segment centers, producing discrete relations

$$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$3

For a shared metasurface interface, the same current unknowns appear in the field operators of both adjoining regions, and the GSTC block supplies the inter-region coupling (Stewart et al., 2018).

This version of LD-BEM is especially well matched to open-region scattering problems because the radiation condition is satisfied naturally through the Green’s functions and no volumetric meshing or absorbing layer is required. It is also the clearest example that LD-BEM does not necessarily imply sparse global matrices: the underlying BEM matrices remain dense, with dense storage and matrix-vector cost $$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$4, although the unknowns are confined to interfaces.

The reported validation problems are both at $$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$5 GHz. In the uniform-metasurface example, the geometry uses a metasurface of length $$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$6 m, a source boundary of length $$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$7 m, a source-to-sheet distance $$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$8 m, and a discretization $$\underline{D}_T^\top \underline{\widetilde{S}}_T \underline{D}_T$$9. LD-BEM results for reflected and transmitted Gaussian-beam fields converge to a semi-analytical plane-wave decomposition. In the nonuniform-metasurface example, a Lorentzian dispersive susceptibility with a modulated resonance frequency $\mathcal{L}$0, using $\mathcal{L}$1 and $\mathcal{L}$2, reproduces diffraction orders whose strengths converge with mesh refinement and agree with an FDFD solution.

6. Numerical characteristics, limitations, and open directions

The computational advantages of LD-BEM depend on which variant is considered. In the time-fractional Fisher–KPP formulation, sparsity is a primary outcome of local-domain partitioning: the assembled global system is sparse and non-symmetric, sparsity increases with the number of subregions, assembly cost scales approximately linearly with the number of subregions, and the memory footprint is dominated by sparse storage of the global matrix together with the fractional history vectors. In the convection-adapted skeletal formulation, sparsity also follows from element-local DtN condensation and skeleton assembly. In the metasurface scattering formulation, by contrast, the global algebra remains dense because the Helmholtz BEM operators are dense, and scalability relies on iterative solvers, preconditioning, and, for large-scale 3D problems, fast multipole or $\mathcal{L}$3-matrix acceleration (Gortsas, 22 Aug 2025, Hofreither et al., 2015, Stewart et al., 2018).

Several limitations recur across the literature. In local-domain diffusion–reaction problems, sparsity and conditioning depend on the number and topology of subregions: overly large subregions reduce sparsity, while excessively small ones increase global degrees of freedom. Strong nonlinearity may exceed the comfortable range of lagged Picard iterations and may require damping or Newton-type updates. Fractional-memory terms introduce a history convolution whose explicit accumulation scales linearly with the number of time steps; memory compression or sum-of-exponentials approximations are identified as plausible improvements but are not implemented in the cited work. In convection-dominated polyhedral problems, very high Peclet numbers may require finer auxiliary face triangulations and more careful face adaptation, and complex polygonal faces make boundary-layer resolution more difficult. In metasurface scattering, EFIE-type first-kind operators are ill-conditioned, especially at low frequency or with dense discretization, so balanced block scaling, Calderón preconditioning, or CFIE-type strategies are relevant.

A second common misconception is that LD-BEM is synonymous with a single physical application area. The documented uses range from anomalous diffusion with logistic reactions, to convection-adapted PDE-harmonic discretizations, to zero-thickness electromagnetic sheets. A plausible implication is that LD-BEM should be classified by its reduction mechanism—local boundary integral or local DtN elimination of interiors—rather than by any single governing equation.

The open directions stated in the cited works are correspondingly broad. For time-fractional diffusion–reaction, natural next steps include higher dimensions, adaptive local domains, advanced convolution quadratures, and integration with hierarchical or FMM accelerations for boundary blocks. For convection-diffusion-reaction, further development is suggested by the absence of new formal a priori error bounds for the convection-adapted traces and by the need to handle more difficult boundary conditions and coefficient variability. For metasurfaces, extensions to fully anisotropic and bi-anisotropic sheets, large-scale 3D RWG/BC discretizations, and improved conditioning machinery are all explicitly compatible with the same GSTC-BEM framework.

In aggregate, LD-BEM occupies a distinctive position among boundary-centered numerical methods. It preserves the core BEM advantage of reducing dimensionality, but it relaxes the strict global boundary-only ideal when nonlinearities, memory terms, convection-dominated layers, or interface constitutive laws make purely classical formulations inefficient or incomplete. Its essential strategy is localize, eliminate, and couple.