# Update on WP2: BEM-HDG coupling

Janosch Preuss

► To cite this version:

Janosch Preuss. Update on WP2: BEM-HDG coupling. Workshop project Butterfly, Sep 2025, Bordeaux, France. hal-05443777

HAL Id: hal-05443777

https://hal.science/hal-05443777v1

Submitted on 6 Jan 2026

HAL is a multi-disciplinary open access
archive for the deposit and dissemination of sci-
entific research documents, whether they are pub-
lished or not. The documents may come from
teaching and research institutions in France or
abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est
destinée au dépôt et à la diffusion de documents
scientifiques de niveau recherche, publiés ou non,
émanant des établissements d'enseignement et de
recherche français ou étrangers, des laboratoires
publics ou privés.

<figure><img src="image_1.png" /></figure>

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
Project Butterfly
Update on WP2: BEM-HDG coupling

Janosch Preuss

Workshop Inria Bordeaux
September 8-10, 2025
Basic idea of the work package
# The solar cycle

* Active regions (e.g. clustering of sunspots)

* Indicate high magnetic activity

* Sun's magnetic field goes through 11 year cycle

* Butterfly diagram

<figure><img src="image_2.png" /></figure>

The data for the synoptic maps is courtesy of the GONG project.
* Solar background model is radially/axi-symmetric ~ can solve efficiently for Green's function

H. Barucq, F. Faucher, D. Fournier, L. Gizon & H. Pham, "Efficient and Accurate Algorithm for the Full Modal Green's Kernel of the Scalar Wave Equation in Helioseismology", SIAP, 2020

* Perturbation $\delta c$ caused by active regions breaks symmetry: How to compute Green's function?

1. Full 3D simulation → WP1

2. **HDG-BEM coupling** → WP2

<figure><img src="image_3.png" /></figure>

Simulations performed by Florian using `hawen` (HDG discretization)
<figure><img src="image_4.png" /></figure>

## Exploit remaining symmetry to reduce computational costs

* **Core assumption:** heterogeneous perturbations $\delta c$ from symmetric background are compactly supported.

* **Main idea:** Solve only in (small) compactly supported regions with expensive 3D method, use Green's function of symmetric case for remainder.

* **Realization:** Coupling of HDG with boundary element method.
* Basic idea of the work package and problem statement (✔)

* Reduction of the problem using boundary integral equations

* Comparison of two HDG-BEM couplings for Helmholtz transmission problem

* Singular quadrature: How to evaluate integrals containing the Green's function?

* Next steps
Reduction of the problem to the scatterers + boundary
using boundary integral formulation
<figure><img src="image_5.png" /></figure>

■ Let $\mathcal{O}$ be the comp. domain for helioseismology (= a ball).

■ Assume active regions contained in $\Omega_{\varepsilon} \in \mathcal{O}$ whereas $\Omega_0 = \mathcal{O} \setminus \bar{\Omega}_{\varepsilon}$.

■ We want to solve for the Green's function

$$(-\Delta - V)G = \delta_s \text{ in } \mathcal{O}, \quad \text{r.b.c. on } \partial\mathcal{O}, \qquad V = V_0 + V_{\varepsilon}$$

■ $V_0$ potential of the (radially symmetric) background model and $V_{\varepsilon}$ represents the perturbation inside the active regions, i.e. $\text{supp } V_{\varepsilon} \in \Omega_{\varepsilon}$.

■ Let $G_0$ denote the Green's function of the background model

$$(-\Delta - V_0)G_0 = \delta_s \text{ in } \mathcal{O}, \quad \text{r.b.c. on } \partial\mathcal{O}.$$
■ We will assume that the source $s$ is located in $\Omega_0$.

■ Introduce new unknowns $(\delta G_0, G_\epsilon)$ such that $G = \begin{cases} G_0 + \delta G_0 & \text{in } \Omega_0, \\ G_\epsilon & \text{in } \Omega_\epsilon. \end{cases}$

■ Now we have to solve

$$ \begin{aligned} (-\Delta - V_0)\delta G_0 &= 0 \text{ in } \Omega_0, &&+ \text{r.b.c. on } \partial\Omega \\ (-\Delta - V)G_\epsilon &= 0 \text{ in } \Omega_\epsilon. \end{aligned} $$

■ With transmission conditions on interface $\Gamma = \bar{\Omega}_0 \cap \bar{\Omega}_\epsilon$

$$ \begin{cases} G_0 + \delta G_0 & = G_\epsilon \text{ on } \Gamma, \\ \partial_n G_0 + \partial_n \delta G_0 & = \partial_n G_\epsilon \text{ on } \Gamma. \end{cases} $$
$$(-\Delta - V_0)\delta G_0 = 0 \text{ in } \Omega_0, \quad \begin{cases} G_0 + \delta G_0 = G_\varepsilon \text{ on } \Gamma, \\ \partial_n G_0 + \partial_n \delta G_0 = \partial_n G_\varepsilon \text{ on } \Gamma. \end{cases} \\ (-\Delta - V)G_\varepsilon = 0 \text{ in } \Omega_\varepsilon.$$

<figure><img src="image_6.png" /></figure>

$$\bullet \mathcal{O} \mapsto \mathbb{R}^d,$$

$$\bullet \text{ r.b.c. on } \partial \mathcal{O} \mapsto \text{Sommerfeld},$$

$$\bullet V_0 \mapsto \omega^2 \text{ constant (!)}$$

$$\bullet \gamma_0 = G_0$$

$$\bullet \Omega_\varepsilon \mapsto \Omega,$$

$$\bullet G_\varepsilon \mapsto u,$$

$$\bullet V \mapsto \omega^2(x),$$

$$\bullet G_0 \mapsto G,$$

$$\bullet \Omega_0 \mapsto \mathbb{R}^d \setminus \Omega,$$

$$\bullet \delta G_0 \mapsto v,$$

$$\bullet \Gamma \mapsto \partial\Omega,$$

$$\bullet \gamma_1 = \partial_n G_0.$$

$$\lim_{|x| \to \infty} \left( \frac{\partial v(x)}{\partial |x|} - i\omega v(x) \right) = 0, \quad \begin{cases} u - v = \gamma_0 \text{ on } \Gamma, \\ \partial_n u - \partial_n v = \gamma_1 \text{ on } \Gamma. \end{cases} \\ (-\Delta - \omega^2(x))u = 0 \text{ in } \Omega.$$

<figure><img src="image_7.png" /></figure>

For Helmholtz in $d=2$ we have $G(x, y) = \frac{i}{4} H_0^1(\omega \|x - y\|)$.
■ There are many ways to do this giving rise to a zoo of different BEM formulations.

► We consider two of them here: symmetric and non-symmetric.

■ Symmetric approach starts from Green's representation formula:

■ Writing $\alpha = v|_{\Gamma}$ and $\mu = \partial_{n}v|_{\Gamma}$ we have for $\tilde{x} \in \mathbb{R}^d \setminus \Omega$

$$ v(\tilde{x}) = - \underbrace{\int_{\Gamma} G(\tilde{x}, y) \mu(y) \, dS(y)}_{:=V\mu(\tilde{x})} + \underbrace{\int_{\Gamma} \partial_{n(y)} G(\tilde{x}, y) \alpha(y) \, dS(y)}_{:=K\alpha(\tilde{x})} $$

■ We need two equations to determine the two unknowns $\alpha$ and $\mu$.

■ Letting $\tilde{x} \to x \in \Gamma$ we obtain the first Calderon identity:

$$ \alpha = -V\mu + \tilde{K}\alpha $$

■ To close the system we use the second Calderon identity:

$$ \tilde{K}^t \mu + W\alpha = 0. $$
■ Single-layer operator

$$V(\varphi)(x) = \int_{\Gamma} G(x, y)\varphi(y) dS(y), \quad x \in \Gamma,$$

■ Double-layer operator

$$K(\varphi)(x) = \int_{\Gamma} \frac{\partial G(x, y)}{\partial \mathbf{n}(y)} \varphi(y) \, dS(y), \quad x \in \Gamma,$$

■ And to ease notation we define $\tilde{K} = \frac{1}{2}I + K$.

■ Hypersingular operator:

$$ \langle W\varphi, \phi \rangle_{\Gamma} = \int_{\Gamma} \int_{\Gamma} \nabla_{\Gamma}\varphi(x)\nabla_{\Gamma}\phi(y)G(x, y)dS(x) \, dS(y) \\ - \omega^2 \int_{\Gamma} \int_{\Gamma} \varphi(x)\phi(y)G(x, y)\mathbf{n}(x) \cdot \mathbf{n}(y)dS(x) \, dS(y) $$
Recall that $\alpha = v|_{\Gamma}$ and $\mu = \partial_{n}v|_{\Gamma}$.

Calderon identities:

$$ \begin{cases} \alpha = -V\mu + \tilde{K}\alpha & \text{on } \Gamma, \\ \tilde{K}^t\mu + W\alpha = 0 & \text{on } \Gamma. \end{cases} $$

Transmission conditions

$$ \begin{cases} u - \alpha = \gamma_0 & \text{on } \Gamma, \\ \partial_n u - \mu = \gamma_1 & \text{on } \Gamma. \end{cases} $$

Using the first Calderon identity to express the continuity condition we obtain

$$ \begin{cases} \partial_n u - \mu = \gamma_1 & \text{on } \Gamma, \\ u + V\mu - \tilde{K}\alpha = \gamma_0 & \text{on } \Gamma, \\ \tilde{K}^t\mu + W\alpha = 0 & \text{on } \Gamma. \end{cases} $$
■ Express exterior field in terms of a density $\rho$:

$$v(\tilde{x}) = K\rho(\tilde{x}), \quad \tilde{x} \in \mathbb{R}^d \setminus \Omega$$

■ Identities for density:

$$\begin{cases} v = \tilde{K}\rho & \text{on } \Gamma, \\ \partial_n v = -W\rho & \text{on } \Gamma. \end{cases}$$

■ Transmission conditions

$$\begin{cases} u - v = \gamma_0 & \text{on } \Gamma, \\ \partial_n u - \partial_n v = \gamma_1 & \text{on } \Gamma. \end{cases}$$

■ Then replacing $v$ and $\partial_n v$ by their density representations in the transmission conditions yields the two conditions:

$$\begin{cases} u - \tilde{K}\rho = \gamma_0 & \text{on } \Gamma, \\ \partial_n u + W\rho = \gamma_1 & \text{on } \Gamma. \end{cases}$$
Symmetric vs. non-symmetric HDG-BEM coupling for
Helmholtz transmission problem
■ To solve the problem in $\Omega$ we use the existing code in `hawen` which is based on a first order formulation of the PDE

$$ \omega \mathbf{q} + \nabla u = 0, \qquad \nabla \cdot \mathbf{q} - \omega u = f, \quad \text{in } \Omega. $$

■ Standard spaces for HDG:

► $V_h^k = \prod_{\tau \in \mathcal{T}_h} [\mathbb{P}_k(\tau)]^d \ni q_h \to q$

► $W_h^k = \prod_{\tau \in \mathcal{T}_h} \mathbb{P}_k(\tau) \ni u_h \to u$

► $M_h^k = \prod_{e \in \mathcal{E}_h} \mathbb{P}_k(e) \ni \lambda_h \to u|_{\mathcal{E}_h}$

■ Spaces for BEM part on $\Gamma = \cup_{\ell=1}^N \Gamma_\ell$:

► $X_h^k := \prod_\ell \mathbb{P}_k(\Gamma_\ell) \ni \mu_h \to \partial_n v|_\Gamma := \mu$

► $Y_h^k = (\prod_\ell \mathbb{P}_k(\Gamma_\ell)) \cap C(\Gamma) \ni \alpha_h \to v|_\Gamma := \alpha$

<figure><img src="image_8.png" /></figure>
$$
\left\{
\begin{array}{lclcl}
\mathbf{q} \cdot \mathbf{n} + \mu & = & -\tilde{\gamma}_1 & \text{on } \Gamma, & \text{Jump condition for flux} \\
\mathbf{u} + V\mu - \tilde{K}\alpha & = & \gamma_0 & \text{on } \Gamma, & \text{Jump condition for trace} \\
\tilde{K}^t \mu + W\alpha & = & 0 & \text{on } \Gamma, & \text{Calderon identity II}
\end{array}
\right.
\quad \rightsquigarrow \quad \text{eq. for } \lambda_h
$$

Find $(\lambda_h, u_h, \mathbf{q}_h, \mu_h, \alpha_h) \in M_h^k \times W_h^k \times \mathbf{V}_h^k \times \mathcal{X}_h^k \times \mathbf{Y}_h^{k+1}$ such that

$$
\begin{align*}
(\omega \mathbf{q}_h, \mathbf{\Psi})_{\mathcal{T}_h} - (u_h, \nabla \cdot \mathbf{\Psi})_{\mathcal{T}_h} + \langle \lambda_h, \mathbf{\Psi} \cdot \mathbf{n} \rangle_{\partial \mathcal{T}_h} &= 0 && \forall \mathbf{\Psi} \in \mathbf{V}_h^k, \\
(\nabla \cdot \mathbf{q}_h, \phi)_{\mathcal{T}_h} - (\omega u_h, \phi)_{\mathcal{T}_h} + (\tau(u_h - \lambda_h), \phi)_{\partial \mathcal{T}_h} &= (f, \phi)_{\mathcal{T}_h} && \forall \phi \in W_h^k, \\
(\tau(u_h - \lambda_h) + \mathbf{q}_h \cdot \mathbf{n}, \xi)_{\partial \mathcal{T}_h} + (\mu_h + \tau_B(\alpha_h - \lambda_h), \xi)_{\Gamma} &= -(\tilde{\gamma}_1 + \tau_B \gamma_0, \xi)_{\Gamma} && \forall \xi \in M_h^k, \\
\langle \lambda_h, \eta \rangle_{\Gamma} + (V \mu_h, \eta)_{\Gamma} - (\tilde{K} \alpha_h, \eta)_{\Gamma} &= (\gamma_0, \eta)_{\Gamma} && \forall \eta \in \mathcal{X}_h^k, \\
\langle \tilde{K}^t \mu_h, \beta \rangle_{\Gamma} + (W \alpha_h, \beta)_{\Gamma} + \langle \tau_B(\alpha_h - \lambda_h), \beta \rangle_{\Gamma} &= -(\tau_B \gamma_0, \beta)_{\Gamma} && \forall \beta \in \mathbf{Y}_h^{k+1}.
\end{align*}
$$
$$ \begin{cases} u - \tilde{K} \rho = \gamma_0 & \text{on } \Gamma, \text{ Jump condition for trace} \rightsquigarrow \text{eq. for } \lambda_h \\ -\mathbf{q} \cdot \mathbf{n} + W \rho = \tilde{\gamma}_1 & \text{on } \Gamma, \text{ Jump condition for flux} \rightsquigarrow \text{eq. for } \varrho_h. \end{cases} $$

Find $(\lambda_h, u_h, \mathbf{q}_h, \varrho_h) \in M_h^k \times W_h^k \times \mathbf{V}_h^k \times \mathbf{Y}_h^{k+1}$ such that

$$ \begin{align*} (\omega \mathbf{q}_h, \mathbf{\Psi})_{\mathcal{T}_h} - (u_h, \nabla \cdot \mathbf{\Psi})_{\mathcal{T}_h} + \langle \lambda_h, \mathbf{\Psi} \cdot \mathbf{\nu} \rangle_{\partial \mathcal{T}_h} &= 0 && \forall \mathbf{\Psi} \in \mathbf{V}_h^k, \\ -(\mathbf{q}_h, \nabla \phi)_{\mathcal{T}_h} - (\omega u_h, \phi)_{\mathcal{T}_h} + (\mathbf{q}_h \cdot \mathbf{\nu} + \tau(u_h - \lambda_h), \phi)_{\partial \mathcal{T}_h} &= (f, \phi)_{\mathcal{T}_h} && \forall \phi \in W_h^k, \\ -(\mathbf{q}_h \cdot \mathbf{\nu} + \tau(u_h - \lambda_h), \xi)_{\partial \mathcal{T}_h \setminus \Gamma} &= 0 && \forall \xi \in M_h^k, \\ \langle \lambda_h, \xi \rangle_{\Gamma} - \langle \tilde{K} \varrho_h, \xi \rangle_{\Gamma} &= (\gamma_0, \eta)_{\Gamma} && \forall \xi \in M_h^k, \\ -(\mathbf{q}_h \cdot \mathbf{n} + \tau(u_h - \lambda_h), \beta)_{\Gamma} + (W \varrho_h, \beta)_{\Gamma} &= (\tilde{\gamma}_1, \beta)_{\Gamma} && \forall \beta \in \mathbf{Y}_h^{k+1}. \end{align*} $$

<footnote>Fu, Heuer, & Sayas, "A non-symmetric coupling of Boundary Elements with the Hybridizable Discontinuous Galerkin method", CAMWA, 2017</footnote>
■ Jump condition for trace

■ Jump condition for flux

Symmetric:

Find $(\lambda_h, u_h, q_h, \mu_h, \alpha_h) \in M_h^k \times W_h^k \times \mathbf{V}_h^k \times \mathbf{X}_h^k \times \mathbf{Y}_h^{k+1}$ such that

$$ (\tau(u_h - \lambda_h) + \mathbf{q}_h \cdot \mathbf{n}, \xi)_{\partial \mathcal{T}_h} + (\mu_h + \tau_B(\alpha_h - \lambda_h), \xi)_\Gamma = -(\tilde{\gamma}_1 + \tau_B\gamma_0, \xi)_\Gamma \quad \forall \xi \in M_h^k, $$

$$ \langle \lambda_h, \eta \rangle_\Gamma + \langle V\mu_h, \eta \rangle_\Gamma - \langle \tilde{K}\alpha_h, \eta \rangle_\Gamma = (\gamma_0, \eta)_\Gamma \quad \forall \eta \in X_h^k, $$

$$ \langle \tilde{K}^t \mu_h, \beta \rangle_\Gamma + \langle W\alpha_h, \beta \rangle_\Gamma + \langle \tau_B(\alpha_h - \lambda_h), \beta \rangle_\Gamma = -(\tau_B\gamma_0, \beta)_\Gamma \quad \forall \beta \in \mathbf{Y}_h^{k+1}. $$

■ The BEM equations (last two) only depend on $\lambda_h \rightsquigarrow$ Given $\lambda_h$ could solve for $(\mu_h, \alpha_h)$.

Non-symmetric:

Find $(\lambda_h, u_h, q_h, \varrho_h) \in M_h^k \times W_h^k \times \mathbf{V}_h^k \times \mathbf{Y}_h^{k+1}$ such that

$$ \langle \lambda_h, \xi \rangle_\Gamma - \langle \tilde{K} \varrho_h, \xi \rangle_\Gamma = (\gamma_0, \eta)_\Gamma \quad \forall \xi \in M_h^k, $$

$$ -(\mathbf{q}_h \cdot \mathbf{n} + \tau(u_h - \lambda_h), \beta)_\Gamma + \langle W\varrho_h, \beta \rangle_\Gamma = (\tilde{\gamma}_1, \beta)_\Gamma \quad \forall \beta \in \mathbf{Y}_h^{k+1}. $$

■ BEM equation for $\varrho_h$ depends directly on $\mathbf{q}_h \cdot \mathbf{n}$ and $u_h$.
■ In both cases: the discontinuous variables $u_h$ and $q_h$ are element-local and can be eliminated from the discrete system by expressing them through $\lambda_h$.

■ Symmetric case $\rightsquigarrow$ linear system for $(\lambda_h, \mu_h, \alpha_h)$

■ Non-symmetric case: $\rightsquigarrow$ linear system for $(\lambda_h, \varrho_h)$.

■ Linear system for non-symmetric case is smaller but has an uglier sparsity structure.

<figure><img src="image_9.png" /><figcaption>Symmetric</figcaption></figure>

<figure><img src="image_10.png" /></figure>

<figure><img src="image_11.png" /><figcaption>Non-symmetric</figcaption></figure>
Test problem: $d = 2$ and $\Omega = \{x \mid \|x\| < 1\}$

Coupled problem:

$$
\left\{
\begin{array}{l}
\omega \mathbf{q} + \nabla u = 0, \\
-\Delta v - \omega^2 v = 0 \text{ in } \Omega_+, \\
u - v = \gamma_0,
\end{array}
\right.
\qquad
\begin{array}{l}
\nabla \cdot \mathbf{q} - \omega u = 0, \quad \text{in } \Omega, \\
\lim_{|x| \to \infty} \left( \frac{\partial v(x)}{\partial |x|} - i\omega v(x) \right) = 0, \\
\mathbf{q} \cdot \mathbf{n} + \partial_n v = -\tilde{\gamma}_1 \quad \text{on } \Gamma.
\end{array}
$$

Sample solution:

$$
\begin{array}{l}
\bullet u(x) = \exp(i\omega x),
\end{array}
$$

$$
\begin{array}{l}
\bullet v(x) = H_0^1(\omega \|x - \tilde{x}\|), \text{ for} \\
\qquad \tilde{x} = (1/5, -1/5)
\end{array}
$$

<figure><img src="image_12.png" /></figure>
Comparison of error $||u - u_h||_{L^2(\Omega)}$ for $k = 2$

■ Both methods achieve the expected convergence rate $||u - u_h||_{L^2(\Omega)} \lesssim h^{k+1}$.

■ Completely similar behavior w.r.t increase of $\omega$.

<figure><img src="image_13.png" /></figure>
- Solid lines (symmetric HDG-BEM coupling)

- Stars (pure BEM): Solve only BEM equations on $\Gamma$ for $(\mu_h, \alpha_h)$ assuming exterior trace $v$ given.

<figure><img src="image_14.png" /></figure>

- Results diverge as $\omega$ increases since HDG suffers noticeably from pollution. Choose higher polynomial degree for HDG than for BEM.
* nzz represents the number of non-zero entries of the original matrix, whereas nzz inv gives the number of entries in the LU factors.

$N = 384$ boundary elements

<table><thead><tr><th>Method</th><th>ndof</th><th>nzz</th><th>nzz inv</th><th>Fact time [s]</th><th>Memory of Fact. [Gb]</th></tr></thead><tbody><tr><td>hdg</td><td>1.26 · 10<sup>5</sup></td><td>2.2 · 10<sup>6</sup></td><td>1.5 · 10<sup>7</sup></td><td>1.9</td><td>0.4</td></tr><tr><td>non-symmetric</td><td>1.27 · 10<sup>5</sup></td><td>4.9 · 10<sup>6</sup></td><td>2.9 · 10<sup>7</sup></td><td>12.9</td><td>0.7</td></tr><tr><td>symmetric</td><td>1.28 · 10<sup>5</sup></td><td>7.5 · 10<sup>6</sup></td><td>2.9 · 10<sup>7</sup></td><td>9.4</td><td>0.8</td></tr></tbody></table>

$N = 768$ boundary elements

<table><thead><tr><th>Method</th><th>ndof</th><th>nzz</th><th>nzz inv</th><th>Fact time [s]</th><th>Memory of Fact. [Gb]</th></tr></thead><tbody><tr><td>hdg</td><td>4.68 · 10<sup>5</sup></td><td>8.4 · 10<sup>6</sup></td><td>6.7 · 10<sup>7</sup></td><td>10</td><td>1.6</td></tr><tr><td>non-symmetric</td><td>4.70 · 10<sup>5</sup></td><td>1.9 · 10<sup>7</sup></td><td>1.2 · 10<sup>8</sup></td><td>92</td><td>3.0</td></tr><tr><td>symmetric</td><td>4.72 · 10<sup>5</sup></td><td>2.9 · 10<sup>7</sup></td><td>1.2 · 10<sup>8</sup></td><td>61</td><td>3.1</td></tr></tbody></table>
* In this example both coupling approaches yield identical accuracy on the same mesh.

* Non-symmetric coupling has significantly less nzz than symmetric coupling.

* Unfortunately this seems not to yield benefits in solution of the linear system: time for factorization as $\omega \uparrow$ is actually lower for symmetric coupling.

* Are the nzz more favourably distributed for the symmetric coupling ?

* In any case: addition of the BEM block to the system increases factorization time drastically ~ consider iterative solution strategies?
Singular quadrature
$$ 
\begin{aligned}
\mathbb{V}_{i,j} &= \langle \mathcal{V} \varphi_j, \varphi_i \rangle_{\Gamma_h} \\
&= \iint_{\Gamma_h} \iint_{\Gamma_h} \mathcal{G}(x, y) \varphi_j(y) \varphi_i(x) \, dx \, dy \\
&= \sum_{n=1}^{N} \sum_{m=1}^{N} \iint_{\Gamma_n} \iint_{\Gamma_m} \mathcal{G}(x, y) \varphi_j(y) \varphi_i(x) \, dx \, dy
\end{aligned}
$$

<figure><img src="image_15.png" /></figure>

* Corresponding matrix $\mathbb{V}$ will be dense due to nonlocality of integral operator.

* Singular integrals for $\bar{\Gamma}_n \cap \bar{\Gamma}_m \neq \emptyset$ as $|G(x, y)| \to \infty$ for $x \to y$

* Requires implementation of special quadrature rules for singular integrals.

* Double integrals over $\Gamma$ and kernel does not live on reference element $\rightsquigarrow$ kernel needs to be evaluated *many* times (dominates cost of assembly)
■ Let us restrict to the two dimensional case.

■ After transforming to the reference configuration we have to approximate integrals of the form

$$ \int_0^1 \int_0^1 f(\hat{x}, \hat{y}) \, d\hat{x} \, d\hat{y} $$

where f has either

<figure><img src="image_16.png" /></figure>

■ a logarithmic singularity at the diagonal $\hat{x} = \hat{y}$ (identical elements $\Gamma_n = \Gamma_m$)

■ a singularity at $(\hat{x}, \hat{y}) = (0, 0)$ (case of adjacent elements upon applying substitution if necessary)

■ or no singularity in the other cases $\rightsquigarrow$ nothing to do.
$$ \int_{0}^{1} \int_{0}^{1} f(\hat{x}, \hat{y}) \, d\hat{x} \, d\hat{y} \approx \sum_{n} \sum_{m} f(\hat{x}_{n,m}, \hat{y}_{n,m}) \omega_{n,m} $$

Vertex singularity

$N_v = 9, M_v = 12$

<figure><img src="image_17.png" /></figure>

Identical elements

$N_s = 10, M_s = 9, p = 5$

<figure><img src="image_18.png" /></figure>
# kernel evaluations (to get down to $10^{-6}$ accuracy)

- On the diagonal need about 5000 evaluations per pair.

- On lower/upper diagonal need less but still a lot.

- Elsewhere the kernel is smooth ~ can get away with lower-order quadrature.

<figure><img src="image_19.png" /></figure>
Next steps
* Replace $G(x, y) = \frac{i}{4}H_0^1(\omega\|x - y\|)$ by numerically computed $G_h(x, y)$,

* Each Green's function evaluation corresponds to one PDE solve.

Questions:

* What is the accuracy of $G_h$ needed to preserve the convergence rates?

* Could also be theoretically investigated for the symmetric approach.

* Will the number of Green function evaluations prove to be too costly?

<figure><img src="image_20.png" /></figure>
- Consider a background media for which $G(x, y)$ is not analytically known.

- First cook up a toy problem (compare with full space solution computed with `hawen`)

- Then we consider the solar background (2d vs. 2.5d ?).
Thank you for your attention!
# Appendix
Symmetric HDG-BEM coupling for Poisson
$$ egin{cases} 
\kappa^{-1}\mathbf{q} + \nabla u = 0, & \nabla \cdot \mathbf{q} = f, \quad \text{in } \Omega, \\
-\Delta v = 0 \text{ in } \Omega_+, & v = O(1/r) \text{ at infinity} \\
u = v + \gamma_0 & -\mathbf{q} \cdot \mathbf{n} = \partial_n v + \gamma_1 \text{ on } \Gamma. 
\end{cases} 
$$

<figure><img src="image_21.png" /></figure>

- The data is assumed to satisfy $\int_{\Omega} f + \int_{\Gamma} \gamma_1 = 0$.

- Symmetric coupling strategy:

Cockburn, Guzmán & Sayas, "Coupling of Raviart-Thomas and Hybridizable Discontinuous Galerkin Methods with BEM", SIAM J. Numer. Anal., 2012
Find $(\lambda_h, p_h, v_h, \mu_h, \alpha_h) \in M_h^k \times W_h^k \times V_h^k \times X_h^k \times Y_h^{k+1}$ such that

$$ (\kappa^{-1}v_h, \Psi)_{\tau_h} - (p_h, \nabla \cdot \Psi)_{\tau_h} + \langle \lambda_h, \Psi \cdot \nu \rangle_{\partial \mathcal{T}_h} = 0 \qquad \forall \Psi \in V_h^k, $$

$$ (\nabla \cdot v_h, \phi)_{\tau_h} + (\tau(p_h - \lambda_h), \phi)_{\partial \mathcal{T}_h} = (f, \phi)_{\tau_h} \qquad \forall \phi \in W_h^k, $$

$$ (\tau(p_h - \lambda_h) + q_h \cdot \nu, \xi)_{\partial \mathcal{T}_h} + (\mu_h + \tau_B(\alpha_h - \lambda_h), \xi)_{\Gamma} = -(\gamma_1 + \tau_B\gamma_0, \xi)_{\Gamma} \qquad \forall \xi \in M_h^k, $$

$$ \langle \lambda_h, \eta \rangle_{\Gamma} + \langle V\mu_h, \eta \rangle_{\Gamma} - \langle \tilde{K}\alpha_h, \eta \rangle_{\Gamma} = (\gamma_0, \eta)_{\Gamma} \qquad \forall \eta \in X_h^k, $$

$$ \langle \tilde{K}^t \mu_h, \beta \rangle_{\Gamma} + \langle W\alpha_h, \beta \rangle_{\Gamma} + \langle \tau_B(\alpha_h - \lambda_h), \beta \rangle_{\Gamma} = -(\tau_B\gamma_0, \beta)_{\Gamma} \qquad \forall \beta \in Y_h^{k+1}. $$

Steps:

(a) Equation for single layer operator (last time: ✅)

(b) Equation for hypersingular layer operator

(c) Calderon identity: couple (a) and (b)

(d) Full HDG-BEM coupling: (a), (b) and all the rest.
(a) Equation for single layer operator (assume $\alpha = v|_{\Gamma}$ given)

Find $\mu_h \in X_h^k$ : $\langle V\mu_h, \eta \rangle_{\Gamma} = \langle (I - \tilde{K})\alpha, \eta \rangle_{\Gamma}, \quad \forall \eta \in X_h^k,$

$$ \begin{align*} \|\mu - \mu_h\|_{H^{-1/2}(\Gamma)} &\lesssim \inf_{x_h \in X_h^k} \|\mu - x_h\|_{H^{-1/2}(\Gamma)} \\ &\lesssim h^{k+3/2} \end{align*} $$

$$ \rightsquigarrow \|\mu - \mu_h\|_{L^2(\Gamma)} \lesssim h^{k+1} $$

■ Expected convergence rates achieved (up to precision of quadrature). ✓

$$ \|\mu - \mu_h\|_{L^2(\Gamma)} $$

<figure><img src="image_22.png" /></figure>
(b) Equation for hypersingular layer operator

Implementational todo's

- $C^0$-elements on $\Gamma$:

$$Y_h^k = \left( \prod_{\ell} \mathbb{P}_k(\Gamma_\ell) \right) \cap C(\Gamma)$$

- Evaluation of derivatives of shape functions along $\Gamma$ for the hypersingular operator:

$$\langle W\varphi, \phi \rangle_{\Gamma} = \int_{\Gamma} \int_{\Gamma} \nabla_{\Gamma}\varphi(x)\nabla_{\Gamma}\phi(y)G(x, y)dS(x) dS(y)$$

- Different polynomial degrees for dG and cG ($X_h^k$ and $Y_h^{k+1}$)
(b) Equation for hypersingular layer operator (assume $\mu = \partial_n v|_\Gamma$ given)

Find $\alpha_h \in Y_h^k$ such that: $$ \langle (W+S)\alpha_h, \beta \rangle_\Gamma = -\langle \tilde{K}^t\mu, \beta \rangle_\Gamma \quad \forall \beta \in Y_h^k. $$

$$ \text{■} (S\alpha_h, \beta)_\Gamma = (\alpha_h, 1)_\Gamma (\beta, 1)_\Gamma $$

$$ \begin{aligned} \|\alpha - \alpha_h\|_{H^{1/2}(\Gamma)} & \lesssim \inf_{y_h \in Y_h^k} \|\alpha - y_h\|_{H^{1/2}(\Gamma)} \\ & \lesssim h^{k+1/2} \end{aligned} $$

$$ \rightsquigarrow \|\alpha - \alpha_h\|_{L^2(\Gamma)} \lesssim h^{k+1} $$

■ Expected convergence rates achieved ✓

<figure><img src="image_23.png" /></figure>
(c) Calderon identity (couple (a) and (b)): the suboptimal way of doing it

Find $(\mu_h, \alpha_h) \in X_h^k \times Y_h^k$ such that

$$ \langle V\mu_h, \eta \rangle_{\Gamma} - \langle \tilde{K}\alpha_h, \eta \rangle_{\Gamma} = -\langle \alpha, \eta \rangle_{\Gamma} \quad \forall \eta \in X_h^k, $$

$$ \langle \tilde{K}^t \mu_h, \beta \rangle_{\Gamma} + \langle (W + S)\alpha_h, \beta \rangle_{\Gamma} = 0 \quad \forall \beta \in Y_h^k. $$

$$ \|\alpha - \alpha_h\|_{H^{1/2}(\Gamma)} + \|\mu - \mu_h\|_{H^{-1/2}(\Gamma)} $$

$$ \lesssim \inf_{y_h \in Y_h^k} \|\alpha - y_h\|_{H^{1/2}(\Gamma)} + \inf_{x_h \in X_h^k} \|\mu - x_h\|_{H^{-1/2}(\Gamma)} $$

$$ \lesssim h^{k+1/2} + h^{k+3/2} $$

$$ \rightsquigarrow \|\alpha - \alpha_h\|_{L^2(\Gamma)} \lesssim h^{k+1}, \quad \|\mu - \mu_h\|_{L^2(\Gamma)} \lesssim h^k $$

$L^2(\Gamma)$ error for $k=2$

<figure><img src="image_24.png" /></figure>

■ Convergence for the flux is one order suboptimal ×
(c) Calderon identity: replace $Y_h^k$ by $Y_h^{k+1}$! $\rightsquigarrow$ optimal rates

$$ \|\alpha - \alpha_h\|_{H^{1/2}(\Gamma)} + \|\mu - \mu_h\|_{H^{-1/2}(\Gamma)} \lesssim \inf_{y_h \in Y_h^{k+1}} \|\alpha - y_h\|_{H^{1/2}(\Gamma)} + \inf_{x_h \in X_h^k} \|\mu - x_h\|_{H^{-1/2}(\Gamma)} \\ \lesssim h^{k+3/2} + h^{k+3/2} $$

$$ \rightsquigarrow \|\mu - \mu_h\|_{L^2(\Gamma)} \lesssim h^{k+1} \qquad \|\alpha - \alpha_h\|_{L^2(\Gamma)} \lesssim h^{k+2} \quad \vrule \\ \|\mu - \mu_h\|_{L^2(\Gamma)} \qquad \|\alpha - \alpha_h\|_{L^2(\Gamma)} $$

<figure><img src="image_25.png" /></figure>

<table><tr><td><pre><code>-* k = 0   ◼ k = 1   ◊ k = 2</code></pre></td><td><pre><code>  O(h<sup>1</sup>)</code></pre></td><td><pre><code>  O(h<sup>2</sup>)</code></pre></td><td><pre><code>  O(h<sup>3</sup>)</code></pre></td></tr></table>

<table><tr><td><pre><code>-* k = 0   ◼ k = 1   ◊ k = 2</code></pre></td><td><pre><code>  O(h<sup>2</sup>)</code></pre></td><td><pre><code>  O(h<sup>3</sup>)</code></pre></td><td><pre><code>  O(h<sup>4</sup>)</code></pre></td></tr></table>
Find $(\lambda_h, p_h, v_h, \mu_h, \alpha_h) \in M_h^k \times W_h^k \times V_h^k \times X_h^k \times Y_h^{k+1}$ such that

$$ (\kappa^{-1} \mathbf{v}_h, \mathbf{\Psi})_{\mathcal{T}_h} - (p_h, \nabla \cdot \mathbf{\Psi})_{\mathcal{T}_h} + \langle \lambda_h, \mathbf{\Psi} \cdot \mathbf{\nu} \rangle_{\partial \mathcal{T}_h} = 0 \qquad \forall \mathbf{\Psi} \in \mathbf{V}_h^k, $$

$$ (\nabla \cdot \mathbf{v}_h, \phi)_{\mathcal{T}_h} + (\tau(p_h - \lambda_h), \phi)_{\partial \mathcal{T}_h} = (f, \phi)_{\mathcal{T}_h} \qquad \forall \phi \in W_h^k, $$

$$ (\tau(p_h - \lambda_h) + \mathbf{v}_h \cdot \mathbf{\nu}, \xi)_{\partial \mathcal{T}_h} + (\mu_h + \tau_B(\alpha_h - \lambda_h), \xi)_{\Gamma} = -(\gamma_1 + \tau_B\gamma_0, \xi)_{\Gamma} \qquad \forall \xi \in M_h^k, $$

$$ \langle \lambda_h, \eta \rangle_{\Gamma} + \langle V\mu_h, \eta \rangle_{\Gamma} - \langle \tilde{K}\alpha_h, \eta \rangle_{\Gamma} = (\gamma_0, \eta)_{\Gamma} \qquad \forall \eta \in X_h^k, $$

$$ \langle \tilde{K}^t \mu_h, \beta \rangle_{\Gamma} + \langle W\alpha_h, \beta \rangle_{\Gamma} + \langle \tau_B(\alpha_h - \lambda_h), \beta \rangle_{\Gamma} = -(\tau_B\gamma_0, \beta)_{\Gamma} \qquad \forall \beta \in Y_h^{k+1}. $$

**Steps:**

(a) Equation for single layer operator (last time: ✔)

(b) Equation for hypersingular layer operator (✔)

(c) Calderon identity: couple (a) and (b) (✔)

**Now:**

(d) Full HDG-BEM coupling: (a), (b) and all the rest.
The 'standard' HDG matrices (in hawen notation)

Vector of dofs: $\mathbf{U} = (p_h, \mathbf{v}_h)$ and $\mathbf{\Lambda} = (\lambda_h)$

$$ \mathbf{A}_e = \begin{pmatrix} \tau(\phi_i, \phi_j)\partial\kappa_e & (\partial_x\phi_i, \phi_j)\kappa_e & (\partial_z\phi_i, \phi_j)\kappa_e \\ -(\phi_i, \partial_x\phi_j)\kappa_e & (\kappa^{-1}\phi_i, \phi_j)\kappa_e & 0 \\ -(\phi_i, \partial_z\phi_j)\kappa_e & 0 & (\kappa^{-1}\phi_i, \phi_j)\kappa_e \end{pmatrix}, \quad \mathbf{S}_e = \begin{pmatrix} (f, \phi_j)\kappa_e \\ 0 \\ 0 \end{pmatrix} $$

$$ \mathbf{C} = \begin{pmatrix} -\tau(\xi_k, \phi_j)_{f_i} \dots -\tau(\xi_k, \phi_j)_{f_{N_{\text{face}}}} \\ (\xi_k, \phi_j \nu_x)_{f_i} \dots (\xi_k, \phi_j \nu_x)_{f_{N_{\text{face}}}} \\ (\xi_k, \phi_j \nu_z)_{f_i} \dots (\xi_k, \phi_j \nu_z)_{f_{N_{\text{face}}}} \end{pmatrix} \mathbf{B}_e = \begin{pmatrix} \mathbf{B}_e^{(f_i)} \\ \vdots \\ \mathbf{B}_e^{(f_{N_{\text{face}}})} \end{pmatrix} \mathbf{L}_e = \begin{pmatrix} \mathbf{L}_e^{(f_i)} \\ \vdots \\ \mathbf{L}_e^{(f_{N_{\text{face}}})} \end{pmatrix} $$

$$ \mathbf{B}_e^{(f_i)} = (\tau_{f_i}(\xi_k, \phi_j)_{f_i}, (\xi_k, \phi_j \nu_x)_{f_i}, (\xi_k, \phi_j \nu_z)_{f_i}) $$

$$ \mathbf{L}_e^{(f_i')} = -\tau_{f_i}(\xi_k, \xi_j)_{f_i}', \quad \mathbf{L}_e^{(f_i^B)} = -(\tau_{f_i} + \tau_B)(\xi_k, \xi_j)_{f_i}^{f_i^B}. $$
$$ \begin{pmatrix} \mathbf{A} & \mathbf{C} & \mathbf{0} & \mathbf{0} \\ \mathbf{B} & \mathbf{L} & (\mu_h, \xi)_{\Gamma} & \tau_B(\alpha_h, \xi)_{\Gamma} \\ \mathbf{0} & (\lambda_h, \eta)_{\Gamma} & \langle V\mu_h, \eta \rangle_{\Gamma} & -\langle \tilde{K}\alpha_h, \eta \rangle_{\Gamma} \\ \mathbf{0} & -\tau_B(\lambda_h, \beta)_{\Gamma} & \langle \tilde{K}^t\mu_h, \beta \rangle_{\Gamma} & \langle (W + \tau_B I)\alpha_h, \beta \rangle_{\Gamma} \end{pmatrix} \begin{pmatrix} \mathbf{U} \\ \mathbf{\Lambda} \\ \mu \\ \alpha \end{pmatrix} = \begin{pmatrix} \mathbf{S} \\ \mathbf{f}_{\lambda} \\ \mathbf{f}_{\mu} \\ \mathbf{f}_{\alpha} \end{pmatrix} $$

■ Note that **U** does only couple to **Λ** but neither to **μ** nor to **α**.

■ Hence, we can form the Schur complement as usual to eliminate $\mathbf{U} = \mathbf{A}^{-1}(\mathbf{S} - \mathbf{C}\mathbf{\Lambda})$ from the equation for $\mathbf{\Lambda}$.
$$ \begin{pmatrix} L - BA^{-1}C & (\mu_h, \xi)_\Gamma & \tau_B(\alpha_h, \xi)_\Gamma \\ (\lambda_h, \eta)_\Gamma & \langle V\mu_h, \eta \rangle_\Gamma & -\langle \tilde{K}\alpha_h, \eta \rangle_\Gamma \\ -\tau_B(\lambda_h, \beta)_\Gamma & \langle \tilde{K}^t\mu_h, \beta \rangle_\Gamma & ((W + \tau_B I)\alpha_h, \beta)_\Gamma \end{pmatrix} $$

<figure><img src="image_26.png" /></figure>

k = 1 and N = 24 bnd elements

<figure><img src="image_27.png" /></figure>

k = 3 and N = 96 bnd elements

<figure><img src="image_28.png" /></figure>
* $\Gamma$ is the unit circle as usual.

* $u = \exp(x) \sin(y)$

* $v = \log \left( \frac{\|x-a\|}{\|x-b\|} \right)$, for $a = (0.5, 0.6)$ and $b = (0.4, 0.4)$.

<figure><img src="image_29.png" /><figcaption>Exact $u$</figcaption></figure>

<figure><img src="image_30.png" /><figcaption>Computed $p_h$ for $k=3$</figcaption></figure>
Expected rates (for exterior variable $v$) are observed:

$$ \|\mu_h - \nabla v \cdot \mathbf{n}\|_{L^2(\Gamma)} \lesssim h^{k+1}, \quad \|\alpha_h - v\|_{L^2(\Gamma)} \lesssim h^{k+2} $$

$$ \|\mu_h - \nabla v \cdot \mathbf{n}\|_{L^2(\Gamma)} $$

$$ \|\alpha_h - v\|_{L^2(\Gamma)} $$

<figure><img src="image_31.png" /></figure>

<table><tr><td>--- k = 1</td><td>○ k = 2</td><td>● k = 3</td></tr><tr><td>--- O(h<sup>2</sup>)</td><td>.... O(h<sup>3</sup>)</td><td>.... O(h<sup>4</sup>)</td></tr></table>

<figure><img src="image_32.png" /></figure>

<table><tr><td>--- k = 1</td><td>○ k = 2</td><td>● k = 3</td></tr><tr><td>--- O(h<sup>3</sup>)</td><td>.... O(h<sup>4</sup>)</td><td>.... O(h<sup>5</sup>)</td></tr></table>
$k = 2$

<table><thead><tr><th>$N \to$</th><th>12</th><th>24</th><th>48</th><th>96</th><th>192</th><th>384</th><th>768</th></tr></thead><tbody><tr><td>ndof (HDG + BEM)</td><td>2.0 · 10<sup>2</sup></td><td>5.8 · 10<sup>2</sup></td><td>2.1 · 10<sup>3</sup></td><td>8.9 · 10<sup>3</sup></td><td>3.4 · 10<sup>4</sup></td><td>1.3 · 10<sup>5</sup></td><td>4.7 · 10<sup>5</sup></td></tr><tr><td>nzz (HDG + BEM)</td><td>7.6 · 10<sup>3</sup></td><td>2.8 · 10<sup>4</sup></td><td>1.1 · 10<sup>5</sup></td><td>4.8 · 10<sup>5</sup></td><td>1.9 · 10<sup>6</sup></td><td>7.6 · 10<sup>6</sup></td><td>2.9 · 10<sup>7</sup></td></tr><tr><td><div style="text-align: center; line-height: 1.2em;"><span style="color: #008080;">ndof (HDG + BEM)</span><br><span style="color: #008080;">ndof (HDG)</span></div><div style="text-align: center; line-height: 1.2em;"><span style="color: #FF0000;">nzz (HDG + BEM)</span><br><span style="color: #FF0000;">nzz (HDG)</span></div></td><td>1.6</td><td>1.3</td><td>1.2</td><td>1.06</td><td>1.03</td><td>1.02</td><td>1.01</td></tr><tr><td>memory solver (Gb)</td><td>1 · 10<sup>-3</sup></td><td>4 · 10<sup>-3</sup></td><td>1.6 · 10<sup>-2</sup></td><td>5.5 · 10<sup>-2</sup></td><td>2.1 · 10<sup>-1</sup></td><td>7.6 · 10<sup>-1</sup></td><td>3.1 · 10<sup>0</sup></td></tr><tr><td><div style="text-align: center; line-height: 1.2em;"><span style="color: #008080;">factorization (s)</span></div></td><td>2 · 10<sup>-3</sup></td><td>8 · 10<sup>-3</sup></td><td>3.6 · 10<sup>-2</sup></td><td>2.0 · 10<sup>-1</sup></td><td>1.2 · 10<sup>0</sup></td><td>7.1 · 10<sup>0</sup></td><td>4.8 · 10<sup>1</sup></td></tr></tbody></table>

$k = 3$

<table><thead><tr><th>$N \to$</th><th>12</th><th>24</th><th>48</th><th>96</th><th>192</th></tr></thead><tbody><tr><td><div style="text-align: center; line-height: 1.2em;"><span style="color: #008080;">ndof (HDG + BEM)</span><br><span style="color: #008080;">ndof (HDG)</span></div><div style="text-align: center; line-height: 1.2em;"><span style="color: #FF0000;">nzz (HDG + BEM)</span><br><span style="color: #FF0000;">nzz (HDG)</span></div></td><td>1.6</td><td>1.3</td><td>1.2</td><td>1.07</td><td>1.03</td></tr><tr><td></td><td>3.9</td><td>4.0</td><td>3.8</td><td>3.3</td><td>3.3</td></tr></tbody></table>

■ Add very few dofs but completely mess up sparsity pattern.
## Numerical verification of HDG-BEM

* Check convergence of $||u - p_h||_\Omega$ and $||q - v_h||_\Omega$ in the interior.

* Preliminary results show $||u - (-p_h)||_\Omega \lesssim h^{k+1}$ (sign-error to be fixed ...)

<figure><img src="image_33.png" /></figure>

## Further steps:

* To be discussed ...

* ?? Same method with numerically computed Green's function ??
Numerical verification of HDG-BEM

■ Sign error was caused by picking the wrong field from ctx_field%field.

- Instead of picking the pressure field(,3), I took the x-component of the gradient field(:,1) which is in this case -p due to the specific choice of the sample solution.

■ Actually, convergence is fine.

$\|u - p_h\|_{L^2(\Omega)}$

<figure><img src="image_34.png" /></figure>

<table>
  <tr>
    <td>k = 1</td>
    <td>k = 2</td>
    <td>k = 3</td>
  </tr>
  <tr>
    <td>Ο(h<sup>2</sup>)</td>
    <td>Ο(h<sup>3</sup>)</td>
    <td>Ο(h<sup>4</sup>)</td>
  </tr>
</table>
Extension of symmetric coupling to Helmholtz
* $G(x, y) = \frac{i}{4} H_0^1(\omega \|x - y\|)$

* Hypersingular operator

$$ \langle W\varphi, \phi \rangle_{\Gamma} = \int_{\Gamma} \int_{\Gamma} \nabla_{\Gamma}\varphi(x)\nabla_{\Gamma}\phi(y)G(x, y)dS(x) \, dS(y) \\ - \omega^2 \int_{\Gamma} \int_{\Gamma} \varphi(x)\phi(y)G(x, y)\mathbf{n}(x) \cdot \mathbf{n}(y)dS(x) \, dS(y) $$

* Note that single layer operator can be obtained essentially 'for free' during assembly of $W$.

* In terms of assembly cost it does not matter whether coupling strategy uses only 2 or all 4 integral operators.

* Careful with complex conjugates ( dot product in Fortran is anti-linear in first component !).
Find $(\lambda_h, p_h, v_h, \mu_h, \alpha_h) \in M_h^k \times W_h^k \times V_h^k \times X_h^k \times Y_h^{k+1}$ such that

$$ (\kappa^{-1}v_h, \Psi)_{\tau_h} - (p_h, \nabla \cdot \Psi)_{\tau_h} + \langle \lambda_h, \Psi \cdot \nu \rangle_{\partial \tau_h} = 0 \quad \forall \Psi \in V_h^k, $$

$$ (\nabla \cdot v_h, \phi)_{\tau_h} - \omega^2(p_h, \phi)_{\tau_h} + (\tau(p_h - \lambda_h), \phi)_{\partial \tau_h} = (f, \phi)_{\tau_h} \quad \forall \phi \in W_h^k, $$

$$ (\tau(p_h - \lambda_h) + v_h \cdot \nu, \xi)_{\partial \tau_h} + (\mu_h + \tau_B(\alpha_h - \lambda_h), \xi)_{\Gamma} = -(\gamma_1 + \tau_B\gamma_0, \xi)_{\Gamma} \quad \forall \xi \in M_h^k, $$

$$ \langle \lambda_h, \eta \rangle_{\Gamma} + \langle V\mu_h, \eta \rangle_{\Gamma} - \langle \tilde{K}\alpha_h, \eta \rangle_{\Gamma} = (\gamma_0, \eta)_{\Gamma} \quad \forall \eta \in X_h^k, $$

$$ \langle \tilde{K}^t \mu_h, \beta \rangle_{\Gamma} + \langle W\alpha_h, \beta \rangle_{\Gamma} + \langle \tau_B(\alpha_h - \lambda_h), \beta \rangle_{\Gamma} = -(\tau_B\gamma_0, \beta)_{\Gamma} \quad \forall \beta \in Y_h^{k+1}. $$

**Calderon identity:** Find $(\mu_h, \alpha_h) \in X_h^k \times Y_h^{k+1}$ such that

$$ \langle V\mu_h, \eta \rangle_{\Gamma} - \langle \tilde{K}\alpha_h, \eta \rangle_{\Gamma} = -\langle \nu, \eta \rangle_{\Gamma} \quad \forall \eta \in X_h^k, $$

$$ \langle \tilde{K}^t \mu_h, \beta \rangle_{\Gamma} + \langle (W+S)\alpha_h, \beta \rangle_{\Gamma} = \langle Sv, \beta \rangle_{\Gamma} \quad \forall \beta \in Y_h^{k+1}. $$
# Pollution: BEM (Calderon identity) vs coupled HDG-BEM ($k=2$)

Shown with marks (circles, crosses) are results for Calderon identity (no coupling).

<figure><img src="image_35.png" /></figure>

BEM does not suffer from pollution, FEM does ~ divergent results as frequency increases

BEM could provide much more accurate results but is restricted by accuracy of the FEM

Should use much higher polynomial degree for HDG than for BEM (otherwise a lot of wasted computational effort)

Accuracy limit of the quadrature probably never reached for moderate-high frequency.
■ We aim to replace $G$ by a numerically computed Green's function $G_h$.

■ How small does $h$ need to be to ensure that accuracy remains unaffected?

■ To test this replace $G$ by $G_\epsilon(x, y) = G(x, y) + \epsilon$ (same for derivative)

$$ \|u - p_h\|_{L^2(\Omega)} $$

<figure><img src="image_36.png" /></figure>

■ Probably will need $G_h$ one order of magnitude more accurate than desired accuracy of $p_h$.
Single equation coupling (Poisson)
Vector of dofs: $\mathbf{U} = (p_h, v_h)$, $\mathbf{\Lambda} = (\lambda_h)$ and $\mathbf{P} = (\varrho_h)$.

$$ \mathbf{A}_e = \begin{pmatrix} \tau(\phi_i, \phi_j)\partial\mathcal{K}_e & (\phi_i\nu_x, \phi_j)\partial\mathcal{K}_e - (\phi_i, \partial_x\phi_j)\mathcal{K}_e & (\phi_i\nu_z, \phi_j)\partial\mathcal{K}_e - (\phi_i, \partial_z\phi_j)\mathcal{K}_e \\ -(\phi_i, \partial_x\phi_j)\mathcal{K}_e & (\kappa^{-1}\phi_i, \phi_j)\mathcal{K}_e & 0 \\ -(\phi_i, \partial_z\phi_j)\mathcal{K}_e & 0 & (\kappa^{-1}\phi_i, \phi_j)\mathcal{K}_e \end{pmatrix}, $$

$$ \mathcal{C} = \begin{pmatrix} -\tau(\xi_k, \phi_j)_{f_1} \dots -\tau(\xi_k, \phi_j)_{f_{N_{\text{face}}}} \\ (\xi_k, \phi_j \nu_x)_{f_1} \dots (\xi_k, \phi_j \nu_x)_{f_{N_{\text{face}}}} \\ (\xi_k, \phi_j \nu_z)_{f_1} \dots (\xi_k, \phi_j \nu_z)_{f_{N_{\text{face}}}} \end{pmatrix} \mathcal{B}_e = \begin{pmatrix} \mathcal{B}_e^{(f_1)} \\ \vdots \\ \mathcal{B}_e^{(f_{N_{\text{face}}})} \end{pmatrix} \mathcal{L}_e = \begin{pmatrix} \mathcal{L}_e^{(f_1)} \\ \vdots \\ \mathcal{L}_e^{(f_{N_{\text{face}}})} \end{pmatrix}, $$
$$-(v_h \cdot \nu + \tau(p_h - \lambda_h), \xi)_{\partial T_h \setminus \Gamma}$$

$$\mathbb{B}_e^{(f_i')} = (\tau_{f_i}(\xi_k, \phi_j)_{f_i}, (\xi_k, \phi_j \nu_x)_{f_i}, (\xi_k, \phi_j \nu_z)_{f_i}), \quad \mathbb{B}_e^{(f_i^B)} = \mathbf{0}$$

$$\mathbb{L}_e^{(f_i')} = \tau(\xi_k, \xi_j)_{f_i'}, \quad \mathbb{L}_e^{(f_i^B)} = (\xi_k, \xi_j)_{f_i^B}.$$

$$-(v_h \cdot \nu + \tau(p_h - \lambda_h), \beta)_{\Gamma}$$

$$\tilde{\mathbb{B}}_e^{(f_i')} = \mathbf{0}, \quad \tilde{\mathbb{B}}_e^{(f_i^B)} = (\tau_{f_i}(\xi_k, \beta_j^{cg})_{f_i}, (\xi_k, \beta_j^{cg} \nu_x)_{f_i}, (\xi_k, \beta_j^{cg} \nu_z)_{f_i}),$$

$$\tilde{\mathbb{L}}_e^{(f_i')} = \mathbf{0}, \quad \tilde{\mathbb{L}}_e^{(f_i^B)} = (\xi_k, \beta_j^{cg})_{f_i'}.$$
$$ \begin{pmatrix} A & C & 0 \\ -B & L & -\tilde{K} \\ -\tilde{B} & \tilde{L} & W+S \end{pmatrix} \begin{pmatrix} U \\ \Lambda \\ P \end{pmatrix} = \begin{pmatrix} S \\ f_{\lambda} \\ f_{\varrho} \end{pmatrix} $$

■ The block $\tilde{B}$ couples the interior variables to the BEM variable $P$.

■ Therefore, need to eliminate $U = A^{-1}(S - CA)$ also from equation for $P$.

$$ \begin{pmatrix} L + BA^{-1}C & -\tilde{K} \\ \tilde{L} + \tilde{B}A^{-1}C & W+S \end{pmatrix} \begin{pmatrix} \Lambda \\ P \end{pmatrix} = \begin{pmatrix} f_{\lambda} + BA^{-1}S \\ f_{\varrho} + \tilde{B}A^{-1}S \end{pmatrix} $$

more complicated sparsity structure of condensed system.

<figure><img src="image_37.png" /></figure>

$k = 1, N = 24$

<figure><img src="image_38.png" /></figure>

$\Lambda$

$P$

$U$
<figure><img src="image_39.png" /></figure>

■ Expected rate for $k = 1, 2$.

■ For $k = 3$ convergence rate is limited by quadrature error.

■ Probably subdivisions necessary to improve accuracy.
Comparison of couplings (Poisson)
<figure><img src="image_40.png" /></figure>

<table><thead><tr><th>Method</th><th>nzz</th><th>Factorization time [s]</th><th>Memory of Factorization [Gb]</th></tr></thead><tbody><tr><td>hdg</td><td>8.4 · 10<sup>6</sup></td><td>15</td><td>1.6</td></tr><tr><td>unsymmetric</td><td>2.9 · 10<sup>7</sup></td><td>88</td><td>3.1</td></tr><tr><td>symmetric</td><td>1.9 · 10<sup>7</sup></td><td>92</td><td>3.0</td></tr></tbody></table>
* Recall: the unsymmetric method is only guaranteed to convergence for $||\kappa||_{L^\infty(\Omega)} < 4$.

* In paper by Sayas et al. is has been seen for a specified example that the condition number can grow drastically (looks like exponential growth) as $\kappa$ increases.

* For our example this is actually not observed:

* Should consider exactly the example from Sayas et al. (BEM code requires some adaptations as bem will be active on more than one facet of the same element).

<figure><img src="image_41.png" /></figure>
Discussion of next steps
- Implement unsymmetric coupling for the Helmholtz equation (+ extension of comparison).

- Condition number studies (example from literature).

- Numerical experiments with different polynomial degrees for HDG and BEM

- Helmholtz with numerically computed Green's function.

- Application to helioseismology (solar coefficients, ...)
■ Pick Dirichlet trace from the harmonic function $u = \exp(x) \sin(y)$.

$$ \|u_1 - t_h\|_{L^2(\Gamma)} $$

<figure><img src="image_42.png" /></figure>
■ Compute the error in terms of number of degrees of freedom (ndof).

$$ \|u_1 - t_h\|_{L^2(\Gamma)} $$

<figure><img src="image_43.png" /></figure>

■ Expected rate of $O(h^{k+1})$ is obtained.

■ Note that higher order is more efficient (lower errors with same ndofs).

■ Quadrature is accurate up to about $10^{-6}$.

■ Going beyond would be difficult.

■ Similar results are obtained with an independent implementation from scratch.
■ Fix mesh with $N = 96$ elements for $k = 2$ and study the error in terms of the number of quadrature points $N_s = M_s$ for identical elements.

■ Compute in double precision (crosses) and quadruple precision (circles).

<figure><img src="image_44.png" /></figure>

■ A large $p$ can lead to faster convergence but may trigger instabilities when double precision is used.

∴ $p = 5$ appears to be a reasonable compromise in double precision.
Boundary elements

<figure><img src="image_45.png" /></figure>
Green's representation theorem

■ Let $-\Delta_y G(x, y) = \delta_y(x)$ and $-\Delta u(y) = 0$ then for $\tilde{x} \in \Omega$:

$$ u(\tilde{x}) = \underbrace{\int G(\tilde{x}, y) \frac{\partial u(y)}{\partial \mathbf{n}(y)} dS(y)}_{:=V(\frac{\partial u}{\partial n})(\tilde{x})} - \underbrace{\int \frac{\partial G(\tilde{x}, y)}{\partial \mathbf{n}(y)} u(y) dS(y)}_{:=K(u)(\tilde{x})} $$

<figure><img src="image_46.png" /></figure>

$$ \begin{cases} -\Delta u = 0 & \text{in } \Omega \\ u = u_0 & \text{on } \Gamma \end{cases} $$

■ For $d=2$ we have e.g. that $G(x, y) = -\frac{1}{2\pi} \log(\|x - y\|)$.

■ For $\tilde{x} \to x \in \Gamma$

$$ V \left( \frac{\partial u}{\partial \mathbf{n}} \right) (x) = \sigma(x) u_0(x) + K(u_0)(x) $$

with $\sigma = 1/2$ if $\Gamma$ sufficiently smooth.
■ We assume from now on that $\Gamma$ is polygonal.

■ Denote $u_0 = u|_{\Gamma}$ and $u_1 = \frac{\partial u}{\partial n}$ and note that $V: H^{-1/2}(\Gamma) \to H^{1/2}(\Gamma)$ and $K: H^{1/2}(\Gamma) \to H^{1/2}(\Gamma)$.

■ Find $u_1 \in H^{-1/2}(\Gamma)$ such that

$$ \langle Vu_1, \tau \rangle_{\Gamma} = \left\langle \left( \frac{1}{2} + K \right) u_0, \tau \right\rangle_{\Gamma} \quad \forall \tau \in H^{-1/2}(\Gamma). $$

■ Partition $\Gamma$ into $N$ boundary elements:

$$ \Gamma = \Gamma_h := \bigcup_{\ell=1}^{N} \Gamma_{\ell}. $$

■ Let $X_h := \prod_{\ell} \mathbb{P}_k(\Gamma_{\ell})$. Note that $X_h \subset H^{-1/2}(\Gamma)$.

■ Find $t_h \in X_h$ such that

$$ \langle Vt_h, \tau_h \rangle_{\Gamma} = \left\langle \left( \frac{1}{2} + K \right) u_0, \tau_h \right\rangle_{\Gamma} \quad \forall \tau_h \in X_h. $$

<figure><img src="image_47.png" /></figure>
$$ \begin{aligned} \langle V\varphi_j, \varphi_i \rangle_{\Gamma_h} &= \iint_{\Gamma_h} \iint_{\Gamma_h} G(x, y)\varphi_j(y)\varphi_i(x) \\ &= \sum_{n=1}^{N} \sum_{m=1}^{N} \iint_{\Gamma_n} \iint_{\Gamma_m} G(x, y)\varphi_j(y)\varphi_i(x) \end{aligned} $$

<figure><img src="image_48.png" /></figure>

* Corresponding matrices **V** and **K** will be dense due to nonlocality of operators

* Singular integrals for $\bar{\Gamma}_n \cap \bar{\Gamma}_m \neq \emptyset$ as $|G(x, y)| \to \infty$ for $x \to y$

* Requires implementation of special quadrature rules for singular integrals.

* Double integrals over $\Gamma$ and kernel does not live on reference element $\rightsquigarrow$ kernel needs to be evaluated many times (dominates cost of assembly)
■ Let us restrict to the Laplace equation in $2d$.

■ After transforming to the reference configuration we have to approximate integrals of the form

$$ \int_0^1 \int_0^1 f(x, y) \, dx \, dy $$

where $f$ has either

<figure><img src="image_49.png" /></figure>

■ a logarithmic singularity at the diagonal $x = y$ (identical elements $\Gamma_n = \Gamma_m$)

■ a singularity at $(x, y) = (0, 0)$ (case of adjacent elements upon applying substitution if necessary)

■ or no singularity in the other cases $\rightsquigarrow$ nothing to do.
Treating the vertex singularity by a Duffy transformation

$$I = \int_0^1 \int_0^1 f(x, y) \, dx \, dy = \underbrace{\int_0^1 \left( \int_0^x f(x, y) \, dy \right)}_{I_1} \, dx + \underbrace{\int_0^1 \left( \int_x^1 f(x, y) \, dy \right)}_{I_2} \, dx$$

■ Substituting $y = ux$ with $0 \le u \le 1$ in $I_1$ we obtain:

$$I_1 = \int_0^1 \int_0^1 f(x, ux) \times dx \, du.$$

■ Treat $I_2$ similarly, after transforming to $0 < x < y$, $0 < y < 1$.

$$\approx I = \int_0^1 \int_0^1 [f(uv, v) + f(v, uv)] \, v \, du \, dv,$$

where weight $v$ will cancel a $1/v$ singularity.

$$I \approx \sum_{n=1}^{N_v} \sum_{m=1}^{M_v} [f(x_m x_n, x_n) + f(x_n, x_m x_n)] x_n \omega_n \omega_m.$$

<figure><img src="image_50.png" /></figure>
Identical elements: derivation (I)

■ The following approach is advocated in the documentation of XLIFE++.

■ First split integral into one part defined on lower triangle $u = x + y < 1$ and second part on upper triangle $x + y > 1$ on line $x + y = 1 + u$ with $0 < u < 1$.

$$I = \int_0^1 \int_0^1 f(x, y) \, dx \, dy = \underbrace{\int_0^1 \left( \int_0^u f(u - v, v) \, dv \right)}_{l_1} \, du + \underbrace{\int_0^1 \left( \int_u^1 f(1 + u - v, v) \, dv \right)}_{l_2} \, du$$

■ After the substitution $v = (1 + s)\frac{u}{2}$ we obtain

$$l_1 = \int_0^1 \frac{u}{2} \int_{-1}^1 f\left((1-s)\frac{u}{2}, (1+s)\frac{u}{2}\right) \, ds \, du,$$

■ which has a logarithmic singularity at $s = 0$.
■ This singularity can be cancelled based on the substitution $s = t^p$ where $p$ is an odd integer.

Kirkup, Yazdani & Papazafeiropoulos, "Quadrature Rules for Functions with a Mid-Point Logarithmic Singularity in the Boundary Element Method Based on the $x = t^p$ Substitution", Am. J. Comput. Math., 2019

■ Hence we obtain

$$I_1 = \int_0^1 \frac{pu}{2} \int_{-1}^1 f \left( (1-t^p)^{\frac{u}{2}}, (1+t^p)^{\frac{u}{2}} \right) t^{p-1} dt du,$$

■ and similarly

$$I_2 = \int_0^1 \frac{p(1-u)}{2} \int_{-1}^1 f \left( \frac{1+u-t^p(1-u)}{2}, \frac{1+u+t^p(1-u)}{2} \right) t^{p-1} dt du,$$
■ Choosing a quadrature rule for $t$ on $(-1, 1)$ that avoids the midpoint ($N_s$ even for Gauss-Lobatto) we obtain:

<figure><img src="image_51.png" /></figure>

$$ I \approx \sum_{m=1}^{M_s} \frac{px_m^u}{2} \omega_m^u \sum_{n=1}^{N_s} \omega_n^t f \left( \left(1 - (x_n^t)^p\right) \frac{x_m^u}{2}, \left(1 + (x_n^t)^p\right) \frac{x_m^u}{2} \right) (x_n^t)^{p-1} \\ + \sum_{m=1}^{M_s} \frac{p(1-x_m^u)}{2} \omega_m^u \sum_{n=1}^{N_s} \omega_n^t f \left( \frac{1+x_m^u-(x_n^t)^p(1-x_m^u)}{2}, \frac{1+x_m^u+(x_n^t)^p(1-x_m^u)}{2} \right) (x_n^t)^{p-1} $$
Non-symmetric coupling with HDG
Fu, Heuer, & Sayas, "A non-symmetric coupling of Boundary Elements with the Hybridizable Discontinuous Galerkin method", CAMWA, 2017

* Express the transmission conditions on $\Gamma$ in terms of a density $\rho \in H^{1/2}(\Gamma)$:

$$u = v,$$

$$-\mathbf{q} \cdot \mathbf{n} = \partial_{\mathbf{n}} v.$$

$$\rightsquigarrow u = \tilde{K}_{\rho},$$

$$\mathbf{q} \cdot \mathbf{n} = W\rho.$$

Find $\mathbf{q}_h \in \mathcal{V}_h$, $u_h \in W_h$, $\hat{u}_h \in M_h$, $\lambda_h$ and $\rho_h \in Y_h$ such that

$$(\kappa^{-1}\mathbf{q}_h, \mathbf{r})_{\mathcal{T}_h} - (u_h, \nabla \cdot \mathbf{r})_{\mathcal{T}_h} + \langle \hat{u}_h, \mathbf{r} \cdot \mathbf{n} \rangle_{\partial \mathcal{T}_h} = 0 \quad \forall \mathbf{r} \in \mathcal{V}_h,$$

$$-(\mathbf{q}_h, \nabla w)_{\mathcal{T}_h} + (\mathbf{q}_h \cdot \mathbf{n} + \tau(u_h - \hat{u}_h), w)_{\partial \mathcal{T}_h} = (f, w)_{\mathcal{T}_h} \quad \forall w \in W_h,$$

$$(\mathbf{q}_h \cdot \mathbf{n} + \tau(u_h - \hat{u}_h), \mu)_{\partial \mathcal{T}_h \setminus \Gamma} = 0 \quad \forall \mu \in M_h,$$

$$\langle \hat{u}_h, \mu \rangle_{\Gamma} - \langle \tilde{K}_{\rho} \rho_h, \mu \rangle_{\Gamma} = 0 \quad \forall \mu \in M_h,$$

$$-(\mathbf{q}_h \cdot \mathbf{n} + \tau(u_h - \hat{u}_h), \rho)_{\Gamma} + \langle W\rho_h, \rho \rangle_{\Gamma} + (\rho_h, \mathbf{1})_{\Gamma}(\rho, \mathbf{1})_{\Gamma} = 0 \quad \forall \rho \in Y_h.$$
* To guarantee that the system has a unique solution (for $h$ sufficiently small) a relevant assumption is

$$ ||\kappa||_{L^{\infty}(\Omega)} < 4. $$

* Assuming sufficient regularity the following error estimate holds

$$ ||\mathbf{q} - \mathbf{q}_h||_{\Omega} + ||\varrho - \varrho_h||_{H^{1/2}(\Gamma)} \lesssim h^{k+1}(|\mathbf{q}|_{H^{k+1}(\Omega)} + ||u||_{H^{k+1}(\Omega)}) + h_{\Gamma}^{k+3/2}||\varrho||_{X^{k+2}(\Gamma)} $$

* Note: density $\varrho$ may not have sufficient regularity near corners of $\Gamma$ 
  ↛ alternative: direct formulation.
- symmetric couplings are based on using two integral equations (four integral operators)

- non-symmetric couplings require only one integral equation (two integral operators)

- main challenge of non-symmetric couplings: obtaining stability

- requires either smoothness assumptions on Γ

- or stringent assumptions on material parameters

## Literature:

Of & Steinbach, "Is the one-equation coupling of finite and boundary element methods always stable?", ZAMM, 2013

> "It turns out that the non-symmetric boundary element approximation of the Steklov-Poincaré operator requires sufficient energy from the interior problem to result in a stable discretization."
Can we make use of an existing BEM solver?
* Need to be able to couple with hawen:

* Exchange of information about discretization on $\Gamma$ needed (map between dofs in hawen and in BEM code)

* In addition to standard FEM-BEM coupling: get Green kernel $G(x, y)$ computed with hawen into BEM library

* should provide high-order shape functions

* should support 2d and 3d

* (mesh format compatible with that used by hawen)
* Developed by Timo Betcke, Matthew Scroggs, ...

* Essentially three different versions:

* Bem++: C++ core with Python-Interface, legacy code 😂, high-order 😊

* Bempp-rs: latest (in early progress 😂) rewrite, written in Rust with Python wrapper

* Bempp-cl: rewrite using PyOpenCL or Numba for just-in-time computation of computation kernels

* tried to hack non-analytic Green kernel into library but did not succeed yet 😂

* coupling to FENICS via respective python interfaces 😂, but haven has no python interface 😂)

* shape functions only up to order 1 😊

* seems to be 3d only 😊
* Developed by Lucy Weggler, Daniel Seibel, Joachim Schöberl

* AddOn to NGSolve ~ cannot be used stand-alone 😊

* High-order shape functions 😊

* 3d-only 😊 (but implementing 2d seems straight-forward)

* Could probably get non-analytic Green kernel $G(x, y)$ in there without much problems 😊 if computed with NGSolve

Currently: trying to implement something similar within hawen
■ For $\lambda = \partial_h v$ and $\varphi = v$ the following identities holds:

$$ \begin{cases} V\lambda - \tilde{K}\varphi = -\varphi, \\ \tilde{K}^t\lambda + W\varphi = 0. \end{cases} $$

■ The corresponding bilinear form

$$ A((\lambda, \phi), (\eta, \psi)) := \langle V\lambda, \eta \rangle_{\Gamma} - \langle \tilde{K}\varphi, \eta \rangle_{\Gamma} + \langle \tilde{K}^t\lambda, \psi \rangle_{\Gamma} + \langle W\varphi, \psi \rangle_{\Gamma}. $$

is coercive on $H^{-1/2}(\Gamma) \times H_0^{1/2}(\Gamma)$ where $H_0^{1/2}(\Gamma) = \{\varphi \in H^{1/2}(\Gamma) \mid \int_{\Gamma} \varphi = 0\}$:

$$ A((\lambda, \phi), (\lambda, \phi)) = \langle V\lambda, \lambda \rangle_{\Gamma} + \langle W\phi, \phi \rangle_{\Gamma} \gtrsim \|\lambda\|_{-1/2, \Gamma}^2 + \|\phi\|_{1/2, \Gamma}^2. $$
Cockburn, Guzmán & Sayas, "Coupling of Raviart-Thomas and Hybridizable Discontinuous Galerkin Methods with BEM", SIAM J. Numer. Anal., 2012

$$V_h^k = \prod_{\tau \in \mathcal{T}_h} [\mathbb{P}_k(T)]^d, \quad W_h^k = \prod_{\tau \in \mathcal{T}_h} \mathbb{P}_k(T),$$

$$X_h^k := \prod_{\ell} \mathbb{P}_k(\Gamma_\ell) \quad Y_h^k = \left( \prod_{\ell} \mathbb{P}_k(\Gamma_\ell) \right) \cap C(\Gamma) \quad M_h^k = \prod_{e \in \mathcal{E}_h} \mathbb{P}_k(e),$$

Find $\mathbf{q}_h \in \underline{V}_h^k$, $u_h \in W_h^k$, $\hat{u}_h \in M_h^k$, $\lambda_h \in X_h^k$ and $\varphi \in Y_h^{k+1}$ such that

$$(\kappa^{-1}\mathbf{q}_h, \mathbf{r})_{\mathcal{T}_h} - (u_h, \nabla \cdot \mathbf{r})_{\mathcal{T}_h} + \langle \hat{u}_h, \mathbf{r} \cdot \mathbf{n} \rangle_{\partial \mathcal{T}_h} = 0 \quad \forall \mathbf{r} \in \mathbf{V}_h^k,$$

$$(\nabla \cdot \mathbf{q}_h, w)_{\mathcal{T}_h} + \langle \tau(u_h - \hat{u}_h), w \rangle_{\partial \mathcal{T}_h} = (f, w)_{\mathcal{T}_h} \quad \forall w \in W_h^k,$$

$$\langle \tau(u_h - \hat{u}_h) + \mathbf{q}_h \cdot \mathbf{n}, \mu \rangle_{\partial \mathcal{T}_h} + \langle \lambda_h + \tau_B(\varphi_h - \hat{u}_h), \mu \rangle_{\Gamma} = 0 \quad \forall \mu \in M_h^k,$$

$$\langle \hat{u}_h, \eta \rangle_{\Gamma} + \langle V\lambda_h, \eta \rangle_{\Gamma} - \langle \tilde{K}\varphi_h, \eta \rangle_{\Gamma} = 0 \quad \forall \eta \in X_h^k,$$

$$\langle \tilde{K}^t \lambda_h, \psi \rangle_{\Gamma} + \langle W\varphi_h, \psi \rangle_{\Gamma} + \langle \tau_B(\varphi_h - \hat{u}_h, \psi) \rangle_{\Gamma} = 0 \quad \forall \psi \in Y_h^{k+1}.$$
■ On $\Gamma$ we have $u = v$

■ and for $\lambda = \partial_n v$ and $\varphi = v$ the Caledron identities

$$ \left\{ \begin{array}{l} V\lambda - \tilde{K}\varphi = -\varphi, \\ \tilde{K}^t\lambda + W\varphi = 0. \end{array} \right. $$

Find $\mathbf{q}_h \in \underline{V}_h$, $u_h \in W_h^k$, $\hat{u}_h \in M_h^k$, $\lambda_h \in X_h^k$ and $\varphi \in Y_h^{k+1}$ such that

$$ (\kappa^{-1}\mathbf{q}_h, \mathbf{r})_{\mathcal{T}_h} - (u_h, \nabla \cdot \mathbf{r})_{\mathcal{T}_h} + \langle \hat{u}_h, \mathbf{r} \cdot \mathbf{n} \rangle \partial_{\mathcal{T}_h} = 0 \quad \forall \mathbf{r} \in \mathbf{V}_h^k, $$

$$ (\nabla \cdot \mathbf{q}_h, w)_{\mathcal{T}_h} + \langle \tau(u_h - \hat{u}_h), w \rangle \partial_{\mathcal{T}_h} = (f, w)_{\mathcal{T}_h} \quad \forall w \in W_h^k, $$

$$ \langle \tau(u_h - \hat{u}_h) + \mathbf{q}_h \cdot \mathbf{n}, \mu \rangle_{\partial \mathcal{T}_h} + \langle \lambda_h + \tau_B(\varphi_h - \hat{u}_h), \mu \rangle_{\Gamma} = 0 \quad \forall \mu \in M_h^k, $$

$$ \langle \hat{u}_h, \eta \rangle_{\Gamma} + \langle V\lambda_h, \eta \rangle_{\Gamma} - \langle \tilde{K}\varphi_h, \eta \rangle_{\Gamma} = 0 \quad \forall \eta \in X_h^k, $$

$$ \langle \tilde{K}^t \lambda_h, \psi \rangle_{\Gamma} + \langle W\varphi_h, \psi \rangle_{\Gamma} + \langle \tau_B(\varphi_h - \hat{u}_h, \psi) \rangle_{\Gamma} = 0 \quad \forall \psi \in Y_h^{k+1}. $$
<caption>TABLE 4.1</caption>

Error estimates for RT-BEM and HDG-BEM in the case of highest regularity. All the estimates assume that $0 < \tau \lesssim 1$ and that $0 \le \tau_B \lesssim h_\Gamma^{-1}$. We have denoted $k^* = k + \min\{k-1, 0\}$.

<table><thead><tr><th></th><th></th><th>RT</th><th>HDG</th><th>Hypotheses:<br>(A) h ≈ h<sub>min,Γ</sub><br>(B) τ<sub>B</sub> ≲ h<sub>Γ</sub> or τ<sub>B</sub> ≈ h<sub>Γ</sub><sup>-1</sup><br>(C) τ ≳ h if k = 0</th></tr></thead><tbody><tr><td>q - q<sub>h</sub></td><td>(L<sup>2</sup>(Ω))<sup>d</sup></td><td>k + 1</td><td>k + 1</td><td></td></tr><tr><td>u - u<sub>h</sub></td><td>L<sup>2</sup>(Ω)</td><td>k + 1</td><td>k + 1</td><td></td></tr><tr><td>ε<sup>u</sup><sub>h</sub></td><td>L<sup>2</sup>(Ω)</td><td>k + 2</td><td>k<sup>*</sup> + 2</td><td></td></tr><tr><td>λ - λ<sub>h</sub></td><td>H<sup>-1/2</sup>(Γ)<br>H<sup>-1</sup>(Γ)<br>L<sup>2</sup>(Γ)</td><td>k + 1<br>k + 3/2<br>k + 1/2</td><td>k + 1<br>k<sup>*</sup> + 3/2<br>k + 1/2</td><td>(A), (B), and (C)<br>(A)</td></tr><tr><td>φ - φ<sub>h</sub></td><td>H<sup>1/2</sup>(Γ)<br>L<sup>2</sup>(Γ)</td><td>k + 1<br>k + 3/2</td><td>k + 1<br>k<sup>*</sup> + 3/2</td><td>(A), (B), and (C)</td></tr><tr><td>ε<sup>u</sup><sub>h</sub></td><td>L<sup>2</sup>(∂T<sub>h</sub>)</td><td>k + 2</td><td>k<sup>*</sup> + 2</td><td>(C)</td></tr></tbody></table>

■ Proven in case $-\Delta v + v = 0$.

■ Γ polyhedral / polygonal

■ $q \in [H^{k+1}(\Omega)]^d$

■ $u \in H^{k+1}(\Omega)$

■ $X^s(\Gamma) = \prod_{\ell=1}^N H^s(\Gamma_\ell)$

■ $\lambda \in X^{k+1/2+s}(\Gamma)$

■ $\varphi \in X^{k+3/2+s}(\Gamma)$, $s \in [0, 1/2]$.