Adaptive Boundary Element Procedure

• Adaptive boundary element procedure is a numerical method that refines the discretization mesh using error indicators to efficiently target singularities and discontinuities.
• It employs a SOLVE–ESTIMATE–MARK–REFINE loop, integrating Galerkin solutions, residual estimators, and Dörfler marking to ensure quasi-optimal convergence rates.
• The approach avoids the saturation assumption by using localized inverse-type inequalities and computable a posteriori error estimators to achieve reliable error control.

An adaptive boundary element procedure is a computational methodology for solving boundary integral equations (BIEs) with rigorously controlled error and computational complexity. Rooted in the framework of the boundary element method (BEM), adaptivity aims to refine the discretization mesh and/or function space selectively, concentrating computational resources in regions where error indicators predict inefficient approximation—typically near geometric singularities, material interfaces, or solution features such as layers or discontinuities. Adaptive strategies for BEM have been formulated for a wide variety of underlying operators, including positive, zero, and negative order (e.g., the Laplace single-layer, double-layer, hypersingular, and combined-field operators), and have been rigorously shown to guarantee quasi-optimal convergence rates without requirement for saturation assumptions. The procedure is algorithmically realized through a nested loop consisting of Galerkin solution, residual-type error estimation, Dörfler marking, and localized mesh/subspace refinement, underpinned by a posteriori error analyses and localized inverse-type inequalities.

1. Abstract Boundary Element Setting and Discretization

The general setting considers BIEs of the form

$Au = f$

where $A:\tilde{H}^t(\Gamma)\to H^{-t}(\Gamma)$ is a linear homeomorphism (self-adjoint and coercive), $\Gamma$ is a polyhedral (possibly open) surface in $\mathbb{R}^3$, $t\in\mathbb{R}$, $u\in\tilde{H}^t(\Gamma)$, $f\in H^{-t}(\Gamma)$. The coercivity and boundedness are characterized by constants $\alpha, \beta>0$: $$\langle Av,v \rangle \geq \alpha \|v\|^2_{\tilde{H}^t(\Gamma)}, \quad \|Av\|_{H^{-t}(\Gamma)} \leq \beta \|v\|_{\tilde{H}^t(\Gamma)}$$ The Galerkin discretization is based on a conforming, shape-regular triangulation $\mathcal{P}$ of $\Gamma$, with piecewise polynomial finite element spaces $S^d_P$ (piecewise constants for $d=0$ or continuous linear for $d=1$) and seeks $u_P\in S^d_P$ such that

$$\langle A u_P, v \rangle = \langle f, v \rangle, \quad \forall v\in S^d_P$$

Galerkin orthogonality and best approximation follow, and the energy norm $|v| = \langle Av, v \rangle^{1/2}$ is equivalent to the Sobolev norm.

2. Adaptive Algorithm: The SOLVE–ESTIMATE–MARK–REFINE Loop

The adaptive BEM is structured around the iterative loop:

1. SOLVE: Compute the Galerkin solution $u_k$ on the current mesh $\mathcal{P}_k$.
2. ESTIMATE: Compute local error indicators $\{\eta_T\}$, constructed from the (possibly oscillation-enriched) residual, and compute the global error estimator $\eta(\mathcal{P}_k)$.
3. MARK: Select a minimal subset $\mathcal{R}_k \subset \mathcal{P}_k$ such that

$$\sum_{T\in\mathcal{R}_k} \eta_T^2 \geq \theta \sum_{T\in\mathcal{P}_k} \eta_T^2$$

for $0<\theta\leq1$ (Dörfler marking).

1. REFINE: Refine marked elements and minimally complete to conformity (e.g., newest-vertex bisection in 2D/3D), ensuring that the global mesh regularity (K-mesh property) is maintained and the complexity bound

$$|\mathcal{P}_{k+1}| - |\mathcal{P}_0| \leq C_c |\mathcal{R}_k|$$

is satisfied.

This loop generates a sequence of nested, admissible triangulations, with the estimator driving local resolution of singularities and error concentrations.

3. A Posteriori Error Estimation and Oscillation

Central to adaptivity is the construction of robust a posteriori error estimators. A unified estimator formula is

$$\eta(\mathcal{P}) = \Big( \sum_{T\in\mathcal{P}} h_T^{2r} \|R_P\|^2_{H^r(T)} \Big)^{1/2}$$

with the residual $R_P = f - A u_P$; $h_T = \mathrm{diam}(T)$; exponent $r$ chosen according to operator order $t$. Local contributions are $\eta_T = h_T^r \|R_P\|_{H^r(T)}$. For improved estimator sharpness and localization, an oscillation term,

$$\mathrm{osc}_r(\mathcal{P}, \omega) = \Big( \sum_{T \subset \omega} h_T^{2r} |R_P|^2_{H^r(T)} \Big)^{1/2}$$

is added, quantifying the portion of the residual unresolved by local polynomial approximation. Estimator reliability (upper bound) and efficiency (lower bound, up to oscillation) are established: $$|u - u_P|^2 + \mathrm{osc}^2(\mathcal{P}, \Gamma) \leq C_1 \eta^2(\mathcal{P})$$

$$\eta^2(\mathcal{P}) \leq C_2 \left( |u - u_P|^2 + \mathrm{osc}^2(\mathcal{P}, \Gamma) \right)$$

for constants $C_1$, $C_2$ independent of mesh or solution.

4. Localized Discrete Bounds and Inverse-Type Inequalities

A critical ingredient is a set of localized, discrete reliability and efficiency estimates that bypass the need for saturation assumptions, essential for optimal adaptivity proofs. For any mesh refinement $\mathcal{P} \preceq \mathcal{P}'$ and the set of refined elements $\Gamma_\star$,

$$|u_{\mathcal{P}} - u_{\mathcal{P}'})|^2 \lesssim \sum_{T\subset \Gamma_\star} \eta_T^2(\mathcal{P}) + \mathrm{osc}^2(\mathcal{P}, \Gamma_\star)$$

$$\sum_{T\subset \Gamma_\star} \eta_T^2(\mathcal{P}) \lesssim |u_{\mathcal{P}} - u_{\mathcal{P}'})|^2 + \mathrm{osc}^2(\mathcal{P}, \Gamma_\star)$$

These bounds hinge on proving a key inverse-type inequality for boundary operators $A$, namely

$$\sum_{T\in\mathcal{P}} h_T^{2(s+t)}\|A v\|^2_{H^s(T)} \leq C_A \|v\|^2_{\tilde{H}^t(\Gamma)}, \quad \forall v\in S_P^d$$

for $s\geq0$ and $s+t$ below the order-regularity threshold. This result is pivotal for localizing error reduction and estimator contraction under refinement.

5. Quasi-Optimal Convergence Rates and Approximation Classes

The adaptive boundary element method achieves geometric error reduction and quasi-optimal algebraic convergence dictated by the best-approximation class,

$$\mathcal{A}_{r,s} = \left\{ u : \forall \varepsilon>0, \; \exists \mathcal{P} \;\mathrm{with}\; |\mathcal{P}|-|\mathcal{P}_0| \lesssim \varepsilon^{-1/s},\; E(\mathcal{P}) \leq \varepsilon \right\}$$

with total error $$E(\mathcal{P})^2 = |u-u_P|^2 + \mathrm{osc}^2(\mathcal{P}, \Gamma)$$. The main convergence theorem states: For any $u\in\mathcal{A}_{r,s}$, the adaptive algorithm (for $0<\theta<\theta^*$): $$E(\mathcal{P}_k) \lesssim (|\mathcal{P}_k|-|\mathcal{P}_0|)^{-s}$$ This establishes that the adaptive method recovers the optimal algebraic rate for $E(\mathcal{P})$ as a function of degrees of freedom, matching the nonlinear approximation class. This result is achieved under shape regularity, local mesh grading control (K-mesh property), and without invoking any saturation assumption on the Galerkin error.

6. Avoidance of the Saturation Assumption

Classical convergence proofs for adaptive BEM and FEM often hinged on the uncomputable "saturation assumption": $$\|u-u_P\| \leq \beta\|u_{P'}-u_P\|, \quad (P' \;\text{uniform refinement of}\; P)$$ The methodology established by Gantumur replaces this with fully computable, localized error/estimator bounds as outlined above. Relying solely on residuals, local oscillations, and elementary mesh data, the analysis mirrors the robust adaptive FEM theory of Dörfler and Morin–Nochetto–Siebert, but extended to the singularly nonlocal setting of BEM.

7. Implementation and Practical Aspects

Algorithmic realization of adaptive BEM with convergence rates involves:

• Conforming and shape-regular mesh management supporting local bisection and refinement, with data structures for both elements and mesh hierarchy.
• Rigorous computation of local error indicators through $H^r$-seminorms, which may require higher-order quadrature, local projection, or difference quotients depending on operator order $t$.
• Oscillation terms, built by polynomial approximation of the local residual, ensure estimator robustness even in the presence of irregular data.
• Marking and refinement routines adhering to Dörfler criteria and minimal closure/completion to maintain conformity.
• Complexity, in the number of degrees of freedom per iteration, scales at most linearly with the number of marks.

Modern implementations leverage efficient storage (e.g., H-matrices for full matrices), localized quadrature, and parallelization at the solve and estimate stages. These methods have been validated for a broad class of boundary integral operators, and cover singularly-perturbed, transmission, and elasticity problems, thus forming a foundational paradigm for high-fidelity, adaptive solution of BIEs in computational and applied mathematics (Gantumur, 2011).