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

HAL
open science

# Gap-SBM: A New Conceptualization of the Shifted Boundary Method with Optimal Convergence for the Neumann and Dirichlet Problems

J Haydel Collins, Kangan Li, Alexei Lozinski, Guglielmo Scovazzi

► To cite this version:

J Haydel Collins, Kangan Li, Alexei Lozinski, Guglielmo Scovazzi. Gap-SBM: A New Conceptualization of the Shifted Boundary Method with Optimal Convergence for the Neumann and Dirichlet Problems. 2025. hal-05206595

HAL Id: hal-05206595

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

Preprint submitted on 11 Aug 2025

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.

Public Domain
# Gap-SBM: A New Conceptualization of the Shifted Boundary Method with Optimal Convergence for the Neumann and Dirichlet Problems

J. Haydel Collins<sup>a</sup>, Kangan Li<sup>c</sup>, Alexei Lozinski<sup>b</sup>, Guglielmo Scovazzi<sup>a</sup>

<sup>a</sup>Department of Civil and Environmental Engineering, Duke University, Durham, North Carolina 27708, USA

<sup>b</sup>UniversitØ Marie et Louis Pasteur, CNRS, LmB (UMR 6623), F-25000 Besançon, France

<sup>c</sup>Energy and Mineral Engineering Department, The Pennsylvania State University, University Park, PA, 16802

## Abstract

We propose and mathematically analyze a new Shifted Boundary Method for the treatment of Dirichlet and Neumann boundary conditions, with provable optimal accuracy in the L²- and H¹-norms of the error. The proposed method is built on three stages. First, the distance map between the SBM surrogate boundary and the true boundary is used to construct an approximation to the geometry of the gap between the two. Then, the representations of the numerical solution and test functions are extended from the surrogate domain to the such gap. Finally, approximate quadrature formulas and specific shift operators are applied to integrate a variational formulation that also involves the elements extended in the gap. An extensive set of two-dimensional tests demonstrates the theoretical findings and the overall optimal performance of the proposed method.

**Keywords:** Shifted Boundary Method; Immersed Boundary Method; small cut-cell problem; approximate domain boundaries; Neumann boundary conditions; untted nite element methods.

## 1. Introduction

This article explores a reimagined approach for the treatment of Dirichlet and Neumann boundary conditions with the Shifted Boundary Method. Although this approach applies to both, the major impetus for this topic stemmed from the pursuit of optimal convergence rates for Neumann conditions in a primal formulation. What follows is a general overview of motivations, historical background, and key ideas that will be explored in depth throughout the present article.

Proper enforcement of Neumann boundary conditions is critical across a wide range of engineering and applied science applications, particularly in computational solid mechanics and heat transfer, where traction and ux boundary conditions are ubiquitous. The analyses performed in these areas rely on computational representations of increasingly complex geometry to the extent that traditional conforming Finite Element Methods may no longer be suiciently practical. These limitations are further exacerbated in situations where meshing must be performed iteratively, as in geometrically demanding scenarios like shape optimization, digital twins, and additive manufacturing.

Addressing the demand for less labor-intensive design algorithms, recent advancements have been made in immersed (or embedded, or untted) computational methods. In the last two decades, much attention has been given on reducing the design cost for problems involving complex geometrical features, described in standard formats (i.e., CAD) and non-standard formats (i.e., STL, level sets, etc.). Indeed, immersed/embedded/untted methods have shown the potential to drastically reduce the pre-processing time involved in the acquisition of the geometry and the generation of the computational grid. An (incomplete) list of these developments in the context of nite element methods include the Immersed Boundary

Email addresses: jhc63@duke.edu (J. Haydel Collins), kbl5610@psu.edu (Kangan Li), alexei.lozinski@univ-fcomte.fr (Alexei Lozinski), guglielmo.scovazzi@duke.edu (Guglielmo Scovazzi)
Finite Element Method (IB-FEM) [9, 61], the cutFEM [6, 16 20, 22 24, 39 42, 50, 53, 55 57, 60], the Finite
Cell Method [36, 40, 41, 54], Immerso-Geometric Analysis [42], B-spline immersed methods [55, 56] and
similar earlier methods. Many of these approaches require the geometric construction of the partial ele-
ments cut by the embedded boundary (cut-cells) to form the solution space, and typically employ Nitsche's
method for consistent weak boundary enforcement.

CutFEM relies on data structures that are considerably more tedious to implement with respect to
corresponding tted nite element methods. Furthermore, integrating the variational forms on the char-
acteristically irregular cut cells may also be di cult and advanced quadrature formulas might need to be
employed [36, 54]. Additionally, small-cut cells can induce poor matrix conditioning and even numerical
instabilities, which need to be addressed with appropriate stabilization operators [15, 21] or element aggre-
gation techniques [6].

The Shifted Boundary Method was proposed as an alternative un tted method that remedies the so called
small cut cell problem by removing cut cells entirely. Instead, the location where boundary conditions
are applied is shifted from the true to an approximate (surrogate) boundary composed of facets belonging
to fully intact cells (the surrogate domain). This shift can be thought of as a modi cation/correction
of the boundary conditions by way of Taylor expansions in the direction of the closest distance between
corresponding points on the true and surrogate boundaries. These shifted conditions are enforced weakly,
using Nitsche's method, leading to a relatively simple, robust, accurate, and e cient algorithm. Indeed, the
computational infrastructure based on distances instead of cut cells can be less algorithmically burdensome
in regards to implementation.

The Shifted Boundary Method (SBM) was introduced in [51] and belongs to the more speci c class of
approximate domain methods [7, 8, 10, 13, 14, 25 27, 38, 49], along with -FEM [30 35] albeit with some
key di erences. The work in Main et al. [51] demonstrated the viability of the SBM for Poisson and Stokes
ow problems. Soon after, the method was generalized in [52] to the advection-di usion and Navier-Stokes
equations, and later to hyperbolic conservation laws in [58]. An analysis of the stability and accuracy of the
SBM for the Poisson, advection-di usion, and Stokes operators was also included in [2, 51, 52],respectively.
A high-order version of the SBM was proposed in [5], applications to solid and fracture mechanics problems
were presented in [3, 44, 46 48] and simulations of static and moving interfaces were developed in [29, 43].
Most recently, the SBM was extended to contact problems in solid mechanics in Li et al. [45]

Until now, Dirichlet and Neumann conditions were handled di erently due to the certain challenges
posed by the limitations of the Taylor expansions. That is to say, the solution and its gradient are available
within piecewise-linear interpolation spaces (P1), however, the higher-order terms in the Taylor expansion
are not. This means that Neumann boundary conditions, if shifted naively, will result in a loss of one order of
convergence in L². Earlier work in Atallah et al. [4] addressed this challenge via a mixed shifted formulation
strategy (e.g. solving for strains and displacements in solid mechanics). Although mixed formulations
increase computational cost, it was shown that it is only necessary to solve the mixed variation along
the strip of elements adjacent to the Neumann boundary while enforcing continuity of the normal stress
component between the primal and mixed domains. This approach successfully recovered second order
accuracy for Neumann Boundary conditions with only a meager increase in computational burden.

Naturally, the potential for a SBM that accurately enforces Neumann boundary conditions without a
mixed formulation remained an open research question. This current work demonstrates the viability of a
newly conceptualized SBM that optimally enforces both Neumann and Dirichlet boundary conditions, which
requires no mixed formulation strategy, and retains the classic Taylor expansion shift operator paradigm.

The key idea in the proposed SBM variant, named here Gap-SBM, is to construct a geometric approx-
imation to the gap between the surrogate and true boundaries, and then devise approximate quadrature
formulas to integrate a modi ed variational formulation, which includes the extension of the solution and
test functions from the surrogate domain to the approximate true boundary.

Although the integration of the variational form in the gap between surrogate and true boundaries
is a prerogative of cutFEM approaches, the proposed method remains conceptually an SBM, because the
construction of the approximation to the gap geometry is done via distance vectors/maps, and the numerical
integration is performed on the surrogate boundary using special quadratures that do not involve any cut
cells. Because these geometric constructions and integration formulas do not involve cut cells, the method
<figure><img src="image_2.png" /><figcaption>Figure 1: The surrogate domain, its boundary, and the distance vector d.</figcaption></figure>

is inherently of SBM type. Specifically, the proposed method needs only the finite element infrastructure that already exists on the original background mesh: no additional degrees of freedom, cut-cells, or ghost penalization are needed.

In the following derivations, analysis, and numerical experiments, the proposed method is shown to be stable and have optimal error convergence rates (in the $L^2$- and $H^1$-norm of the error). We emphasize that, to the best of our knowledge, it is for the first time in the literature that optimal convergence in the $L^2$ norm can be proven theoretically for a method acting on untied grids without using advanced quadrature formulas on cut cells. Indeed, previous analyses of different SBM and -FEM variants had to contend with a half-order sub-optimality, which was never observed in practice, but was persistent theoretically. Here, the fully optimal convergence is proven, at least for the symmetric variant of the Gap-SBM.

The rest of this article is organized as follows: Section 2 introduces the SBM notation, Section 3 derives the new proposed SBM for the Poisson problem with Dirichlet and Neumann conditions, Section 4 derives stability results and error estimates in the $H^1$- and $L^2$-norms, Section 5 extends the proposed method to the equations of compressible isotropic linear elasticity, and Section 6 demonstrates the optimal convergence of the error in a series of numerical tests.

## 2. Preliminaries on the Shifted Boundary Method

This section introduces the notation and general strategy of the Shifted Boundary Method (SBM). In Section 3, we will describe the specific details of the new SBM pursued in this work.

### 2.1. Surrogate domains and boundaries

Let be a connected open set in $R^2$ with Lipschitz boundary = @ and let n be the outer-pointing normal to D. We consider a closed domain D such that $\text{clos}(D)$ and we introduce a family $T_h$ of admissible and shape-regular, quasi-uniform tessellations (i.e., grids, or meshes) of D. We will indicate by $h_T$ the size of element $T$ and by h the piecewise constant function such that $h_{jT} = h_T$. In the numerical experiments, we will consider tessellations that are either triangular or Cartesian. For triangular grids, shape-regularity is intended in the sense of Ciarlet. For Cartesian grids, D is uniformly discretized with square elements of side h. The numerical analysis will be restricted to triangular grids to avoid the complex notation and proofs for the general setting. In this context, the quasi-uniformity hypothesis is reframed, with a slight abuse of notation, setting the function h to be globally constant.
Remark 1. The assumption of quasi-uniformity is not essential for the numerical analysis of the proposed methods, but it greatly simplifies the notation in the mathematical proofs.

Remark 2. In this work we limit the discussion to two dimensions, but analogous strategies can be applied in the three-dimensional case. We leave this extension to future, more applied work.

As shown in Figure 1, the SBM is based on restricting the tessellation where the discrete variational formulation is applied to those elements that are strictly contained in $clos()$, i.e., we form

$$T_h := fT \ 2 \ T_h : T \ clos() \ g;$$

which identifies the surrogate domain

$$\begin{array}{r@{\quad}l}
0 & 1 \\
\sim_h := \text{int} \ @ \ [ & T A \\
& T 2 T_h
\end{array};$$

with surrogate boundary $\sim_h := @\sim_h$ and outward-oriented unit normal vector $\mathbf{n}$ to $\sim_h$. Obviously, $T_h$ is an admissible and shape-regular tessellation of $\sim_h$ (see Figure 1a). In other words, all cut elements are removed from the active computational domain, which is now $\sim_h$ instead of $\emptyset$. We now introduce a mapping

$$M_h : \sim_h ! \quad ; \qquad (1a)$$

$$\mathbf{x} \ 7! \ \mathbf{x} \quad ; \qquad (1b)$$

which associates to any point $\mathbf{x} \ 2 \ \sim_h$ on the surrogate boundary a point $\mathbf{x} = M_h(\mathbf{x})$ on the physical boundary. Whenever uniquely defined, the closest-point projection of $\mathbf{x}$ upon $\mathbf{x}$ is a natural choice for $\mathbf{x}$, as shown e.g. in Figure 1b. Through $M_h$, a distance vector function $d_{M_h}$ can be defined as

$$d_{M_h}(\mathbf{x}) = \mathbf{x} \quad \mathbf{x} = [M_h \quad I](\mathbf{x}) : \qquad (2)$$

For the sake of simplicity, we set $d = d_{M_h}$ where $d = \text{kdk}$ and $\mathbf{x}$ is a unit vector.

Remark 3. If $\mathbf{x} = M_h(\mathbf{x})$ does not belong to corners or edges, then the closest-point projection implies $= \mathbf{n}$, where $\mathbf{n}$ has been defined as the outward pointing normal to $\emptyset$.

Remark 4. There are strategies for the definition of the map $M_h$ and distance $d$ other than the closest-point projection, such as level sets, for which $d$ is defined by means of a distance function. Other more sophisticated choices of $M_h$ may be locally preferable and we refer to [1] for more details.

In case the boundary is partitioned into a Dirichlet boundary $D$ and a Neumann boundary $N$ with $= \overline{D} \ [ \ N \text{ and } D \ \setminus \ N = ; \text{, we need to identify whether a surrogate edge } \mathbf{e} \ \sim_h \text{ is associated with } D \text{ or } N \text{. To that end, we partition } \sim_h \text{ as } \overline{\sim_h;D} \ [ \ \overline{\sim_h;N} \text{ with } \overline{\sim_h;D} \ \setminus \ \overline{\sim_h;N} = ; \text{ using again a map } M_h, \text{ such that}$

$$\sim_{h;D} = f\mathbf{e} \quad \sim_h : M_h(\mathbf{e}) \quad D g \qquad (3)$$

and

$$\sim_{h;N} = f\mathbf{e} \quad \sim_h : M_h(\mathbf{e}) \quad N g : \qquad (4)$$

We will also assume that $\sim_{h;N} = \sim_h n \ \overline{\sim_{h;D}}$, that is, that either a surrogate edge $\mathbf{e}$ entirely belongs to the surrogate Dirichlet boundary or to the surrogate Neumann boundary. We will then prevent the case of mixed Dirichlet/Neumann surrogate edges. This hypothesis is realized in practice, by renouncing to the map $M_h$ as the closest-point projection, as discussed in more detail in [1].
## 2.2. General notation for inner products, norms, and seminorms

Throughout this article, we denote by $L^2(!)$, for $! \ge 0$, the space of Lebesgue square-integrable functions on $[0, T]$. We will use the Sobolev spaces $H^m(!) = W^{m;2}(!)$ of index of regularity $m \ge 0$ and index of summability 2, equipped with the (scaled) norm

$$kv_{H^m(!)} = kv k_{L^2(!)}^2 + \sum_{k=1}^{n} \frac{k1(!)^k D^k v k_{L^2(!)}^2}{!^{1-2}}; \quad (5)$$

where $D^k$ is the kth-order spatial derivative operator and $l(A) = \text{meas } 2(A)^{1=2}$ is a characteristic length of the domain A. Note that $H^0(!) = L^2(!)$ and, as usual, we use a simplified notation for norms and semi-norms, i.e., we set $kv_{H^m(!)} = kv_{H^m(!)}$ and $jv_{k;1} = kD^k v_{k_0} = kD^k v_{k_{L^2(!)}}$.

We also introduce the definition of the L²-inner product over $\mathbb{R}$, namely $(u; v)_! = \langle u, v \rangle$, and an analogous inner product on the subset $\subsetneq$ , namely $hu; wi = uw$. We can also restrict to $! \ge 0$ and the norms and seminorms initially defined on $\subsetneq$, that is $k_{k;1} \ge k_{j;1}$ and $k_{k;0} \ge k_{j;0}$, for example.

## 2.3. General strategy of the standard Shifted Boundary Method

Consider now the Poisson problem with Dirichlet and Neumann boundary conditions:

$$u = f; \quad \text{in } \Omega; \quad (6a)$$

$$u = u_D; \quad \text{on } \partial D; \quad (6b)$$

$$r u n = h_N; \quad \text{on } \partial N; \quad (6c)$$

where $\overline{D} = \{ \overline{N} \text{ and } \partial N = \emptyset \}$. Note also that vectors and tensors are marked in bold, while scalars are marked with regular fonts.

As already mentioned, the SBM discretizes the governing equations in $\tilde{h}$ rather than in $h$, with the challenge of accurately imposing boundary conditions on $\tilde{h}$. To this end, boundary conditions are shifted from $\tilde{h}$ to $h$, by performing an mth-order Taylor expansion of the variable of interest at the surrogate boundary, under the assumption that a solution variable $u$ is suitably smooth in the strip between $\tilde{h}$ and $h$. Let $D_d^i$ denote the ith-order directional derivative along $d$:

$$D_d^i u = \sum_{2 \le n \le j} \frac{i!}{j!} \frac{\partial u}{\partial x^j} d;$$

Then, for $x \in \tilde{h}$ and $x = M_h(x)$ and we can write

$$u(x) = u(x+d(x)) = u(x) + \sum_{i=1}^{n} \frac{D_d^i u(x)}{i!} + (R^m(u; d))(x); \quad (7)$$

where the remainder $R^m(u; d)$ satisfies $jR^m(u; d)j = o(\text{kd}^m)$ as $\text{kd} \ge 0$. Assume that the Dirichlet condition $u(x) = u_D(x)$ needs to be imposed on the true boundary $\partial D$. Using the map $M_h$, one can extend $u_D$ from $\partial D$ to $\tilde{h};\partial D$ as $\text{Eu}_D(x) = u_D(M_h(x))$. Then, the Taylor expansion can be used to enforce the Dirichlet condition on $\tilde{h};\partial D$ rather than $\partial D$, as

$$S_d^m u + \text{Eu}_D + R^m(u; d) = 0; \quad \text{on } \tilde{h};\partial D; \quad (8)$$

where we have introduced the boundary shift operator for every $x \in \tilde{h};\partial D$, namely:

$$S_d^m u(x) := u(x) + \sum_{i=1}^{n} \frac{D_d^i u(x)}{i!}; \quad (9)$$
Neglecting the remainder Rm(u; d), we obtain the nal expression of the shifted approximation of order m
of the boundary condition

$$
S_d^m u \quad \text{Eu}_D; \quad \text{on } \sim_{h;D}^\gamma \tag{10}
$$

This shifted boundary condition will be enforced weakly in what follows, and whenever there is no source
of confusion, the symbol E will be removed from the extended quantities, and we would write u<sub>D</sub> in place
of Eu<sub>D</sub>.

$$
S_d u \quad u_D; \text{ on } \sim_{h;D}^\gamma; \tag{11}
$$

with

$$
S_d := S_d^1 u(\star) = u(\star) + r u(\star) d(\star); \quad (12)
$$

To develop a simple SBM variational formulation, assume a triangular grid and consider the space

$$
V_h(\sim_h) = v_h 2 C^0(\sim_h) j v_{hjT} 2 P^1(T); 8T 2 T_h; \quad (13)
$$

where P¹(T) is the the space of linear polynomials over the triangle T 2 T~h. The penalty-free SBM formulation inspired by the work in [28] reads:

Find $u_h 2 V_h(\sim_h)$ such that, $8w_h 2 V_h^k(\sim_h)$

$$
\begin{equation}
\begin{split}
& (r u_h; r w_h)_{\sim_h} \quad \text{hr} \ u_h \ \mathbf{r}; w_h i_{\sim_h} \quad \text{hn} \ \mathbf{r} (h_N \ r \ u_h \ n); w_h i_{\sim_{h:N}} \\
& \qquad + h S_d u_h \quad u_D; r w_h \ \mathbf{r} i_{\sim_{h;D}} \quad (f; w_h)_{\sim_h} = 0:
\end{split}
\tag{14}
\end{equation}
$$

It relies on Nitsche's method to enforce Dirichlet boundary conditions, in which the shift operator $S_d$ (i.e., a Taylor expansion in our case) is used to extend the discrete solution $u_h$ $2 V_h(\sim_h)$ from $\sim_{h;D}$ to $D$. Neumann conditions are imposed weakly assuming a constant extrapolation of the gradient $r u_h$. This choice is due to the fact that, in the case of piecewise-linear approximation, it is not possible to construct a Taylor expansion of $r u_h$, since the Hessian of $u_h$ would vanish. It was found in numerical computations that this limitation causes the SBM to have suboptimal convergence (i.e., $\text{rst-}$ instead of second-order convergence) in the $L^2$-norm of the error, while the $H^1$-seminorm of the error converges optimally (\text{rst-order}). Adapting the discrete approximation spaces, similar derivations and conclusions can be obtained in the case of bi-linear elements over Cartesian grids.

Atallah et al. [4] addressed the sub-optimality of the SBM with Neumann boundary conditions via a mixed shifted formulation (e.g. solving for strains and displacements in solid mechanics). In the next section, we will develop an alternative approach based on a primal formulation.

3. A new conceptualization of the Shifted Boundary Method

To cure the loss of optimality in the convergence of the L²-norm, we propose to account for the e ect
on the solution u_h of the gap n ~h between ~h and (the gray-shaded region in Figure 1a). Of course,
we will avoid explicitly integrating the discrete equations over the gap region but, rather, we will derive
approximate quadrature formulas that do not require integration over cut elements.

To this end, we define the extension $V_h^{\text{ext}}()$ of $V_h(\sim_h)$, where a function $v_h^{\text{ext}}$ $2 V_h^{\text{ext}}()$ is obtained as the linear combination of the extensions to the gap $n \sim_h$ of the piecewise-linear basis functions used to represent $v_h$ in $\sim_h$. In other words, referring to the sketch of Figure 2 for the a two-dimensional triangular grid, each shape function that is non-zero over the elements attached to edges in $\sim_h$ is evaluated over a point in the gap, then the linear combination forming $v_h$ is taken and renamed $v_h^{\text{ext}}$. For all other elements in the discretization we have instead that $v_h^{\text{ext}} = v_h$, and no explicit extension is needed.

The gap $n \sim_h$ is discretized as follows: for every element with an edge $e$ on $\sim_h$ (e.g., the element $T_1$ in Figure 2), we consider the two end nodes of that edge (e.g., $\alpha_1$ and $\alpha_2$), and project them via $M_h$ onto (to $a_1^{\text{ext}} = M_h(\alpha_1)$ and $a_2^{\text{ext}} = M_h(\alpha_2)$), to obtain then a projection of $e$, called $e^{\text{ext}}$, which interpolates between $a_1^{\text{ext}}$ and $a_2^{\text{ext}}$. We define the union of the edges $e^{\text{ext}}$ as $h$, the interpolant of $e$, which can be
<figure><img src="image_3.png" /><figcaption>Figure 2: The extension of the shape functions from the surrogate domain (e.g., elements $T_1$ and $T_2$) to the gap $n \sim_h$ (e.g., the element extensions $T_1^{\text{ext}}$ and $T_2^{\text{ext}}$). $T_1^{\text{ext}}$ has the edge $e^{\text{ext}}$ that interpolates $n \sim_h$, and similarly for $T_2^{\text{ext}}$.</figcaption></figure>

further decomposed into a Dirichlet part $D$; and Neumann part $N$; and $h$. The quadrilateral that connects $a_1$, $a_2$, $a_1^{\text{ext}}$, and $a_2^{\text{ext}}$ is the element extension ($T_1^{\text{ext}}$ in Figure 2) to $n \sim_h$ of the original element ($T_1$). Let us denote by $T_h^{\text{ext}}$ the set of quadrilaterals constructed to discretize the gap $n \sim_h$ by this procedure.

**Remark 5.** In the geometric construction presented here, the boundary is approximated as polygonal, that is $h$, where $h$ is a polygonal domain with its vertices lying on the boundary of $\sim_h$. $h$ introduces a geometric error, which is however quadratic in nature, and this approximation can be made safely in the context of piecewise-linear finite element approximation spaces. Also observe that one can avoid this approximation by using the distance $d$ along the entire edge $e$, so that the edge $e^{\text{ext}}$ would be curved. We prefer to avoid these complications for the sake of simplicity. In the case of higher-order discretizations, the previous argument can be adjusted by computing distances at each node along the edges in $\sim_h$ (including nodes internal to the edges), and constructing with such distances a higher-order approximation $h$ of $\sim_h$.

Based on the discussion in Remark 5, we will always use interchangeably $\sim_h$ and $h$ in what follows. The previous geometric construction allows us to extend the shape functions defined over $T_1$ to $T_1^{\text{ext}}$. This is in a nutshell the construction of a function $v_h^{\text{ext}}$ $2$ $V_h^{\text{ext}}()$. In what follows, we will only consider the two-dimensional setting, but analogous derivations can be extended to the three-dimensional case.

Observe that the dimension of the function spaces $V_h^{\text{ext}}$ ($\sim_h$) and $V_h$ ($\sim_h$) are the same, since no additional degrees of freedom are added in the extension process.

We want now to derive a weak form that implements a SBM discretization of the strong form (6). For a suitably regular solution $u$ of the finite element dimensional problem, the strong form (6) is multiplied by $w_h^{\text{ext}}$ $2$ $V_h^{\text{ext}}()$:

$$ (w_h^{\text{ext}}; u) = (w_h^{\text{ext}}; f); \qquad (15) $$

Let us introduce the following notation, for the sake of brevity:

$$ (\cdot;)_T = \begin{array}{c} X \\ T \cdot \end{array}; \qquad (16a) $$

$$ (\cdot;)_h = \begin{array}{c} X \\ T \cdot \end{array}; \qquad (16b) $$
$$ h; i_{T_h} = \begin{cases} X & h; i_{\mathcal{T}}; \\ \mathcal{T}^2 \mathcal{T}_h & \end{cases} \qquad (16c) $$

$$ h; i_{\mathcal{T}_h^{\text{ext}}} = \begin{cases} \mathcal{T}^2 \mathcal{T}_h & h; i_{\mathcal{T}^{\text{ext}}}; \\ \mathcal{T}^{\text{ext}} 2 \mathcal{T}_h^{\text{ext}} & \end{cases} \qquad (16d) $$

Then integrating by parts equation (15) over all the triangles in $\mathcal{T}_h$ and the quadrilaterals in $\mathcal{T}_h^{\text{ext}}$, we have

$$ \begin{aligned} & (r w_h; r u)_{\mathcal{T}_h} h w_h; r u n_T i_{\mathcal{T}_h} + (r w_h^{\text{ext}}; r u)_{\mathcal{T}_h^{\text{ext}}} h w_h^{\text{ext}}; r u n_{\mathcal{T}^{\text{ext}}} i_{\mathcal{T}_h^{\text{ext}}} \\ & = (w_h; f)_{\mathcal{T}_h} + (w_h^{\text{ext}}; f)_{\mathcal{T}_h^{\text{ext}}}; \end{aligned} \qquad (17) $$

where we have used the fact that for a given $w_h$ and $V_h(\sim_h)$, the extension $w_h^{\text{ext}}$ coincides with $w_h$ over $\sim_h$.

**Remark 6.** While $w_h$ and $V_h(\sim_h)$ are continuous over $\sim_h$, the extension $w_h^{\text{ext}}$ and $V_h^{\text{ext}}$ can be discontinuous over the gap $n\sim_h$. Looking at Figure 2, and in particular to the edge $e_{\alpha_2}$ emanating from $\alpha_2$ along the distance vector $d(\alpha_2)$, it is clear that the extended shape functions from the elements $\mathcal{T}_1$ and $\mathcal{T}_2$ may not match on $e_{\alpha_2}$. This happens because the gradient over $\mathcal{T}_1^{\text{ext}}$ is the same as the one over $\mathcal{T}_1$ and the gradient over $\mathcal{T}_2^{\text{ext}}$ is the same as the one over $\mathcal{T}_2$, but the gradients over $\mathcal{T}_1$ and $\mathcal{T}_2$ in general do not match.

Because of the presence of potential discontinuities, the framework of Discontinuous Galerkin Methods seems the most appropriate to proceed. The contribution from the internal element boundaries in equation (17) can be expanded by making use of the following notations and identities of the jumps and averages of edge quantities:

$$ \|w\| = w^+ n^+ + w n; \qquad (18) $$

$$ \|v\| = v^+ n^+ + v n; \qquad (19) $$

$$ f w \mathbf{g} = \frac{1}{2} (w^+ + w); \qquad (20) $$

$$ \|v w\| = f v \mathbf{g} \|w\| + \|v\| f w \mathbf{g}; \qquad (21) $$

In what follows, to simplify the notation, we will write $w_h$ in place of $w_h^{\text{ext}}$ and $V_h^{\text{ext}}$, since the regions of integration uniquely define whether we are considering a test/trial function on $\sim_h$ or its extension over $n\sim_h$. Hence, denoting $E^o$ the set of interior faces in $\sim_h$ and $E^{\text{ext}, o}$ the set of faces of elements in $\mathcal{T}_h^{\text{ext}}$ that do not lie on nor $\sim_h$, we have

$$ \begin{aligned} & h w_h; r u n i_{\mathcal{T}_h} + h w_h; r u n i_{\mathcal{T}_h^{\text{ext}}} = h 1; \| [w_h r u] \| i_{E^o [ \sim_h [ E^{\text{ext}, o} ]]} \\ & \qquad + h w_h; \int_{\mathcal{Z}} u_n i_{N;h} + h w_h; r u n i_{D;h} \\ & = h 1; f w_h \mathbf{g} \| [r u] \| + \| [w_h] \| \text{ fr } u \mathbf{g} i_{E^o [ \sim_h [ E^{\text{ext}, o} ]]} \\ & \qquad + h w_h; h_N i_{N;h} + h w_h; r u n i_{D;h} \\ & = h \| [w_h] \|; \text{ fr } u \mathbf{g} i_{E^{\text{ext}, o}} + h w_h; h_N i_{N;h} \\ & \qquad + h w_h; r u n i_{D;h}; \end{aligned} \qquad (22) $$

where we have used the fact that $\| [r u] \| = 0$ on edges in $E^o [ \sim_h [ E^{\text{ext}, o} ]$ and $\| [w_h] \| = 0$ on edges in $E^o [ \sim_h]$. Observing that $(r w_h; r u)_{\mathcal{T}_h} = (r w_h; r u)_{\sim_h}$ and $(r w_h; f)_{\mathcal{T}_h} = (r w_h; f)_{\sim_h}$, we have

$$ \begin{aligned} & (r w_h; r u)_{\sim_h} + (r w_h; r u)_{\mathcal{T}_h^{\text{ext}}} h \| [w_h] \|; \text{ fr } u \mathbf{g} i_{E^{\text{ext}, o}} h w_h; h_N i_{N;h} \\ & h w_h; r u n i_{D;h} = (w_h; f)_{\sim_h} + (w_h; f)_{\mathcal{T}_h^{\text{ext}}}; \end{aligned} \qquad (23) $$
Let us now replace the infinite dimensional exact solution $u$ with its approximation $u_h^{\text{ext}}$ and simplify the notation with $u_h$ in place of $u_h^{\text{ext}}$ and $V_h^{\text{ext}}$, as we did before for $w_h^{\text{ext}}$. Complementing the previous equation with terms that weakly enforce the Dirichlet condition and the continuity of the solution across edges in $E^{\text{ext}}$, we obtain the discrete weak formulation:

$$ (r w_h; r u_h)_{\sim h} + (r w_h; r u_h)_{T_h^{\text{ext}}} h w_h; h_N i_{N;h} \\ h [[w_h]]; \text{fr } u_h g i_{E^{\text{ext}}; o} h w_h; r u_h n i_{D;h} \text{ h fr } w_h g; [[u_h]] i_{E^{\text{ext}}; o} \text{ h r } w_h n; u_h u_D i_{D;h} \\ + h h \, {}^1[[w_h]]; [[u_h]] i_{E^{\text{ext}}; o} + h h \, {}^1 w_h; u_h u_D i_{D;h} = (w_h; f)_{\sim h} + (w_h; f)_{T_h^{\text{ext}}}; \quad (24) $$

where Remark 6 is in order also for $u_h$. The terms

$$ h \text{ fr } w_h g; [[u_h]] i_{E^{\text{ext}}; o} \text{ h r } w_h n; u_h u_D i_{D;h} + h h \, {}^1[[w_h]]; [[u_h]] i_{D;h} + h h \, {}^1 w_h; u_h i_{D;h} $$

are typical of an interior penalty discontinuous Galerkin discretization. In particular, for $= 1$ and $> 0$ we obtain a symmetric interior penalty Galerkin discretization, while for $= 1$ and $= 0$ we obtain the skew-symmetric (or non-symmetric) interior penalty Galerkin discretization, which has the advantage of being penalty-free. An SBM version of the latter had been recently explored in [28].

Now, rather than computing integrals on the gap between $\sim h$ and $h$, or on $D;h$ and $N;h$, we propose an approximation to these integrals using information on the $\sim h$ and appropriate rescaling of integrals. Let us start approximating the term:

$$ (r w_h; r u_h)_{T_h^{\text{ext}}} \stackrel{\text{X}}{\sim} h r w_h; r u_h \frac{j T_e^{\text{ext}} j}{j e j} i_e; \quad (25) $$

where, $T_e$ is the quadrilateral emanating from edge $e$. In practice, we have reduced the integrals over extended elements $T^{\text{ext}}$ and $T_h^{\text{ext}}$ to integrals over edges $e \sim h$, by introducing the rescaling factor $j T^{\text{ext}} j = j e j$, which is the ratio between the measure (area) of $T^{\text{ext}}$ and the measure (length) of the edge $e$. See Figure 3 for a sketch of the geometric construction. This simplified approach is reminiscent of using a biased left-node quadrature over an interval in one space dimension, instead of the common mid-point or Gauss quadratures. Similarly, defining

$$ H_e := \frac{j T_e^{\text{ext}} j}{j e j}; \quad (26) $$

we obtain

$$ (w_h; f)_{T_h^{\text{ext}}} \stackrel{\text{X}}{\sim} h w_h; f H_e i_e; \quad (27) $$

Consider now the term

$$ h w_h; h_N i_{N;h} \stackrel{\text{X}}{\sim} h S_d(w_h); h_N (M_h(\star)) j_e i_e; \quad (28) $$

where we approximated the Jacobian of the transformation mapping any boundary edge $e \sim h$ to its extension $e^{\text{ext}}$ as

$$ j_e := \frac{j e^{\text{ext}} j}{j e j} \frac{j d M_h(\star) j}{j d x j}; \quad (29) $$

thus neglecting a higher-order error contributions. Note also that the map $M_h$ is used to evaluate $h_N$ on $x = M_h(\star)$, and that the shift $S_d(w_h)$ is used to evaluate $w_h$ on $e^{\text{ext}}$. Hence, we use the reference edge $e$ to perform the integration over $e^{\text{ext}}$, introducing an approximate formula for the change of variables under
integration. Similarly:

$$h w_h; r u_h n i_{D,h} \qquad X \qquad h S_d(w_h); (r u_h n) j_e i_e; \qquad (30)$$

$$hr w_h n; u u_D i_{D,h} \qquad e^{\sim h:D} X \qquad hr w_h n; (S_d(u_h) u_D) j_e i_e; \qquad (31)$$

$$hh^1 w_h; u_h u_D i_{D,h} \qquad e^{\sim h:D} X \qquad hh^1 S_d(w_h); (S_d(u_h) u_D) j_e i_e; \qquad (32)$$

$e^{\sim h:D}$

Note that in the approximations (30) (31) we introduce a slight abuse of notation by writing the normal $n$ inside the integrals over $e$. In fact, at any point $x \in 2e$, the normal vector $n$ to the actual boundary should be evaluated as $n(M_h(x))$, and the same applies $u_D$ in (31) (32). The last three remaining terms to be approximated are:

$$h[[w_h]]; \text{fr } u_h g_i E^{\text{ext}, o} \qquad X \qquad \begin{array}{l} jd_{a,j} [[S_{1=2d} w_h]]_{\mathbf{k}} \text{ fr } u_h g_{a,j} \\ \text{a} \cdot 2 \text{N } (\sim_h) \\ = X \qquad \frac{jd_{a,j}}{2} [[r w_h d_{a,j}]_{\mathbf{k}} \text{ fr } u_h g_{a,j}]; \end{array} \qquad (33)$$

$$h \text{fr } w_h g; [[u_h]]_{E^{\text{ext}, o}} \qquad X^{\text{a} \cdot 2 \text{N } (\sim_h)} \qquad jd_{a,j} \text{fr } w_h g_{a,j} [[S_{1=2d} u_h]]_{\mathbf{k}} \qquad (34)$$

$$hh^1 [[w_h]]; [[u_h]]_{E^{\text{ext}, o}} \qquad X^{\text{a} \cdot 2 \text{N } (\sim_h)} \qquad = \frac{jd_{a,j}}{2} \text{fr } w_h g_{a,j} [[r u_h d_{a,j}]_{\mathbf{k}}; \qquad (34)$$

$X^{\text{a} \cdot 2 \text{N } (\sim_h)}$

$$h h^1 [[w_h]]; [[u_h]]_{E^{\text{ext}, o}} \qquad X^{\text{a} \cdot 2 \text{N } (\sim_h)} \qquad \begin{array}{l} \bar{h} \, jd_{a,j} [[S_{1=2d} w_h]]_{\mathbf{k}} [[S_{1=2d} u_h]]_{\mathbf{k}} \\ = \frac{jd_{a,j}}{4h} [[r w_h d_{a,j}]_{\mathbf{k}} [[r u_h d_{a,j}]_{\mathbf{k}}; \end{array} \qquad (35)$$

$X^{\text{a} \cdot 2 \text{N } (\sim_h)}$

where $S_{1=2d}(w_h(x)) = w_h(x) + 1 = 2r w_h(x)$ and $N(\sim_h)$ is the set of grid nodes along the surrogate boundary $\sim_h$, at which both $u_h$ and $w_h$ are continuous. In particular, the average $f_{g_a}$ and jump $[[\mathbf{k}]_a$ across the node $\mathbf{a}$, $\sim_h$ are computed by evaluating the elds on the two edges of $\sim_h$ emanating from $\mathbf{a}$, with formulas analogous to (18) and (19), assuming that the unit vectors $n^+$; $n$ are normal to the edge in $E^{\text{ext}, o}$ stemming from $\mathbf{a}$, i.e. to the vector $d_{\mathbf{a}}$. In practice, formulas (33) (35) approximate the integrals on the edges in $E^{\text{ext}, o}$ with a mid-point type formula, using information available at each node $\mathbf{a}$, $2 \cdot N(\sim_h)$. The magnitude $jd_{a,j}$ represents the length of each edge. Note that the formulas (33) (34) are actually exact in the case of triangular meshes and linear finite elements since the gradients of involved functions are then constant on each cell. Approximations (33) (35) can also be applied to the case of Cartesian grids with bi-linear elements, and should include the higher-order terms in the Taylor expansion in the case of higher-order finite element spaces. In summary, the weak form 24 is approximated as

$$(r w_h; r u_h)_{\sim_h} + X \qquad X \qquad h S_d(w_h); h N(M_h(x)) j_e i_e \\ e^{\sim_h} \qquad e^{\sim_h:N} X \qquad \frac{jd_{a,j}}{2} [[r w_h d_{a,j}]_{\mathbf{k}} \text{ fr } u_h g_{a,j} + \text{fr } w_h g_{a,j} [[r u_h d_{a,j}]_{\mathbf{k}}] \\ \text{a} \cdot 2 \text{N } (\sim_h) X \qquad h S_d(w_h); (r u_h n) j_e i_e + \text{hr } w_h n; (S_d(u_h) u_D) j_e i_e \\ e^{\sim_h:D}$$
<figure><img src="image_4.png" /><figcaption>Figure 3: An edge e ∈ E, the attached element T_e ∈ T, and the attached extended element T_e^{ext} ∈ T. The above sketch also features the geometric construction for the integration of the variational forms on the extended element T_e^{ext}. Observe that the positions of the quadrature points on the lateral edges of T_e^{ext} are q_1^{ext} = a_1 + 1 = 2d_{a_1} and q_2^{ext} = a_2 + 1 = 2d_{a_2}, respectively. Instead, q_{e_1} and q_{e_2} are Gauss quadrature points along the edge e.</figcaption></figure>

$$ 
\begin{aligned}
& + \sum_{\substack{\mathbf{x} \\ a=2 \text{ N}(\sim_h)}} \overline{4h} j d_{\mathbf{x}} j [\mathbf{r} w_h d_{\mathbf{x}}]_{\mathbf{k}} [\mathbf{r} u_h d_{\mathbf{x}}]_{\mathbf{k}} + \sum_{\mathbf{e} \sim_h} \overline{h} h S_d(w_h) (S_d(u_h) u_D) j_{\mathbf{e}} i_{\mathbf{e}} \\
& = (w_h; f)_{\sim_h} + \sum_{\mathbf{e} \sim_h} h w_h; f H_{\mathbf{e}} i_{\mathbf{e}} : \quad (36)
\end{aligned}
$$

Remark 7. Observe that the shape functions utilized by the proposed method are a partition of unity. Furthermore, the proposed formulation is exact if the solution is a ne, because the gradient is globally constant and the Taylor expansion is exact in this case.

## 4. Theoretical analysis

As already alluded to at the beginning of Section 2.1, the theoretical analysis will be performed only for two-dimensional triangular grids. With some appropriate adjustment, the analysis could be extended to uniform Cartesian grids, but we prefer to omit this step for the sake of brevity and to avoid complex notation. For the sake of simplicity and without loss of generality, we shall only consider here problem (6) with homogeneous boundary conditions, that is the case of $u_D = 0$ and $h_N = 0$. Problem (6) is well posed for any $f \in L^2(\Omega)$, and we assume that the boundary is suiently smooth to have elliptic regularity, that is the weak solution $u$ is assumed in $H^2(\Omega)$ with $juj_2; 6 Ckf k_0;$. This choice is motivated by the fact that we would ultimately like to study the $L^2$-optimality of the numerical approximation $u_h$ to $u$ in the case of Dirichlet and Neumann boundary conditions, and this can be done only in the case in which $u$ is at least in $H^2(\Omega)$.

As already mentioned in Remark 5, in principle the quadrilaterals $T_e^{ext}$ have curved edges $e^{ext}$, which are approximated with straight edges between nodes. The union of these straight edges forms the boundary $h = @^h$ of the approximate domain $h \sim_h$. The measure of $h$ approximates the measure of up to an error $O(h^2)$. Neglecting the discrepancy between $h$ and is normally acceptable in the case of piecewise-linear (globally continuous) nite element spaces, since the solution error in the natural and $L^2$ norms are typically of order $O(h)$ and $O(h^2)$, respectively. Hence, for the moment, we neglect this ne level of approximation and develop our theory as if one can compute the length of the edge $je^{ext}$ j and the area
of the quadrilateral $jT_h^{\text{ext}}$ j exactly. Moreover, the Jacobian of the mapping from an edge $e$ on the surrogate boundary to the corresponding edge $e^{\text{ext}}$ is approximated according to (29). Taking these errors into account would make the theoretical analysis much more tedious. While this is a necessity in the case of higher-order approximation spaces, these errors can be neglected for the case of piecewise-linear interpolation spaces considered here. In the analysis that follows we will consider the following bilinear form,

$$ 
\begin{aligned}
a_h(w; v) ={}& (r w; r v)_{\sim_h} + (r w; r v)_{T_h^{\text{ext}}} \\
& + h [[w]]_X \text{ fr } v \mathfrak{g} i_{E^{\text{ext}}, o} + h [[v]]_X \text{ fr } w \mathfrak{g} i_{E^{\text{ext}}, o} + h h^{\text{ext}} [[w]]_X [[v]]_{E^{\text{ext}}, o} \\
& + h S_d(w); S_d(r v) n j_e i_e + h S_d(r w) n; S_d(v) j_e i_e \\
& + e^{\sim_h} X^{h; D} + h h^{\text{ext}} S_d(w); S_d(v) j_e i_e; \\
& + e^{\sim_h} X^{h; D}
\end{aligned} 
\quad (37) $$

which is an intermediate step between (24) and (36), in the sense that it relates to (36) but, to simplify the notation, the terms (33), (34), and (35) are not approximated and (29) is applied as for (30), (31), and (32). Here, we interpret the shift operator $S_d$ as

$$ S_d(v^h) = v^h(M_h(x)); \quad (38) $$

an idea we already explored in [28, 59]. When applied to a piecewise linear function $v_h$ 2 $V_h$, this definition of $S_d$ is consistent with (i.e., identical to) definition (12), which is based on Taylor expansions. Moreover, for $v_h$ 2 $V_h$, $S_d(r v_h) = r v_h$ on any edge $e \sim_{h;D}$, since $r v_h$ is piecewise constant. Thus the SBM solution $u_h$ 2 $V_h$ of problem (6) with homogeneous boundary conditions satisfies

$$ a_h(v_h; u_h) = (v_h; f)_{\sim_h} + \begin{aligned}[t] & X \\ & e^{\sim_h} \end{aligned} + \begin{aligned}[t] & \text{h}v_h; \text{f} H_e i_e; & 8v_h \ 2 \ V_h: \end{aligned} \quad (39) $$

Interpreting the shift operator $S_d$ as in (38) allows us to write that the exact solution $u$ 2 $H^2(\Omega)$ to the same problem satisfies

$$ a_h(v_h; u) = (v_h; f); \quad 8v_h \ 2 \ V_h; \quad (40) $$

where we also have neglected errors in the geometric approximation of the boundaries according to (29), so that

$$ 
\begin{aligned}
X \quad & h S_d(v_h); S_d(r u) \quad n j_e i_e = \begin{aligned}[t] & X \quad Z \\ & e^{\sim_{h;D}} \end{aligned} + S_d(v_h) S_d(r u) \quad n \frac{\text{jd}M_h(x)j}{\text{jd}xj} \\
& e^{\sim_{h;D}} X \quad Z \\
& = \begin{aligned}[t] & v_h \quad r \quad u \quad n \\ & e^{\sim_{h;D}} M_h(e) \end{aligned} \\
& = \text{h}v_h; \text{r} \quad u \quad n i_{0;\sim_{h;D}}: \end{aligned} 
\quad (41) $$

An important assumption is now in order:

Assumption 1. The term $j_e$ is bounded above, that is $jj_e j \in C_D$.

Observe that Assumption 1 is normally satisfied in engineering computations, even the ones with the most challenging geometries. This is because the grids utilized in practical computations produce a surrogate boundary that broadly captures the shape of the true domain, for any grid resolution.

The first results in our analysis are two lemmas on the stability of the bilinear form (39), followed by an error estimate in the natural norm.

Lemma 1 (stability of the penalty-based formulations). The bilinear form $a_h$ is coercive on $V_h$ for any 2 $\mathbb{R}$ and $> 0$ sufficiently large. More precisely, there exist $> 0$ and $> 0$ such that
regularity of the mesh and on , such that for all > 0 and all v_h 2 V_h

$$
a_h(v_h; v_h) > 9 v_h 9_h^2
$$

with

$$
9 v_h 9_h^2 := j v_h j_1^2 + \frac{1}{h} k[[v_h]] k_{0;E \text{ ext}; o}^2 + \frac{1}{h} k S_d(v_h) k_{0; \sim_{h;D}^{|j|}}^2 \quad (42)
$$

$$
\begin{align}
& kwk_{0;\sim_{h;D}^{|j|}} := \int_X kwk_{0;j[e]} ; \tag{43} \\
& \phantom{kwk_{0;\sim_{h;D}^{|j|}} := \int_X kwk_{0;j[e]} ; } e \sim_{h;D} \int_X \dots \nonumber
\end{align}
$$

$$
\text{kwk}_{0; j; [e]} := w \frac{p}{j_e} \int_{0; e} : \quad (44)
$$

Proof

We start by stating a trace inverse inequality (which implicitly relies on Assumption 1)

$$
k\mathrm{f}\nabla v_h gk_{0;E\mathrm{ext};o}^2 + kr v_h k_{0;\sim_{h;D}^{|j|}}^2 \geq 6 \frac{C_{\mathrm{inv}}^2}{h} jv_h j_1^2;
$$

and we conclude

$$
\begin{align*}
a_h(v_h; v_h) &> j v_h j_1^2 && (1+\frac{1}{X}) h[[v_h]] ; \text{ fr } v_h g i_{E \text{ ext}; o} + h h \frac{1}{X} [[v_h]] ; [[v_h]] i_{0;E \text{ ext}; o} \\
&\quad + (1+\frac{1}{X}) h S_d(v_h) ; r v_h n j_e i_{0;e} &&+ h h \frac{1}{X} S_d(v_h) ; S_d(v_h) j_e i_{0;e} \\
&\qquad e \sim_{h;D} && e \sim_{h;D} \\
&> 1 \frac{1}{2} C_{\text{inv}}^2 j v_h j_1^2 &&+ \frac{1}{h} \frac{(1+\frac{1}{2})^2}{2''h} k[[v_h]] k_{0;E \text{ ext}; o}^2 + k S_d(v_h) k_{0;\sim_{h;D}^{|j|}}^2 :
\end{align*}
$$

This gives the announced estimate taking " sufficiently small, e.g. " = C<sub>inv</sub><sup>2</sup>, and 0 = <sup>1+(1+<sup>2</sup>)<sup>2</sup>C<sub>inv</sub><sup>2</sup></sup>&frasl;<sub>2</sub>, so that
= <sup>1</sup>&frasl;<sub>2</sub>.

The above lemma does not cover the antysymmetric penalty-free case, i.e. = 1, = 0. Indeed, the bilinear form a_h is coercive with respect to the H^1 seminorm, but not the norm 9 9_h. However, we can replace the notion of coercivity by the inf-sup property, as stated in the following lemma. This can be done under additional assumption about the Dirichlet part of the surrogate boundary:

Assumption 2. There exist positive constants c_D and C_D such that on any edge e ~h;D we have c_D 6
je_6 C_D and n > c_D, where, at any point * 2 e, n is the normal to e that points outside of ~h, while n
is the normal to the actual boundary evaluated at M_h(*).

Observe that condition $\mathbf{n} > c_D$ in Assumption 2 can be violated, but typically only on isolated edges. We conjecture that a numerical analysis could then help to prove the upcoming result even without this assumption.

Lemma 2 (Stability of the penalty-free formulation ). Assume = 1, = 0. There exists > 0,
depending only on the regularity of the mesh, such that

$$
\inf_{v_h \ge 2} \sup_{w_h \ge 2} \frac{a_h(w_h; v_h)}{9 w_h \ 9_h \ 9 v_h \ 9_h} > :
$$

Proof

Taking w_h = v_h, we immediately obtain the coercivity of a_h with respect to the H^1 seminorm:

$$
a_h(v_h; v_h) > j v_h j_1^2 ; \quad 8 v_h 2 v_h :
$$

Now, for any $v_h 2 V_h$, introduce $g_h$ as the piecewise-linear function on $T_h$ taking the same values as $v_h$ at
the boundary nodes on ~h;D and set to zero on at all the other nodes of the mesh. Reexamining the proof of [28, Lemma 1], we observe that, for any boundary edge e ~ h;D,

$$
hr g_h \mathbf{n} ; S_d(v_h)j_e i_e = hr g_h \mathbf{n}(\mathbf{n} \cdot \mathbf{n}) ; S_d(v_h)j_e i_e + hr g_h \sim(\sim \mathbf{n}) ; S_d(v_h)j_e i_e \\
> c_D \frac{c_1}{h} k v_h \frac{p}{j_e} k_{0;e}^2 \quad c_D c_2 C_D j v_h j_{1;T_e}^2 \quad jhr v_h \sim(\sim \mathbf{n}) ; S_d(v_h)j_e i_e j : \quad (46)
$$

where $T_e$ denotes the element of $T_h$ attached to $e$. This bound corresponds to equation (14a) in the proof of [28, Lemma 1], taken in an element-wise version (prior to summing over all the boundary edges). In order to adapt (45) to the needs of the current article, in particular to pass from the normal vector $\mathbf{n}$ to $\mathbf{n}$, we need to invoke Assumption 2. Introducing the unit tangent vector $\mathbf{t}$ on $e$ alongside the normal vector $\mathbf{n}$, we can derive from (45):

$$
\begin{align}
hr g_h \mathbf{n} ; S_d(v_h)j_e i_e ={}& hr g_h \mathbf{n}(\mathbf{n} \cdot \mathbf{n}) ; S_d(v_h)j_e i_e + hr g_h \sim(\sim \mathbf{n}) ; S_d(v_h)j_e i_e \nonumber \\
& > c_D \frac{c_1}{h} k v_h \frac{p}{j_e} k_{0;e}^2 \quad c_D c_2 C_D j v_h j_{1;T_e}^2 \quad jhr v_h \sim(\sim \mathbf{n}) ; S_d(v_h)j_e i_e j : \tag{46}
\end{align}
$$

In the last term above we have replaced $r g_h \sim$ by $r v_h \sim$ since $g_h = v_h$ on $e$. By scaling arguments and Young inequality, we can further bound this term as

$$
jhr v_h \sim(\sim \mathbf{n}) ; S_d(v_h)j_e i_e j 6 \frac{c}{h} \overline{jv_h j_{1;T_e}} k S_d(v_h) k_{0;j} [e] 6 \frac{C}{h} j v_h j_{1;T_e}^2 + c_D \frac{c_1}{4h} k S_d(v_h) k_{0;j}^2 [e] \quad (47)
$$

with h-independent constants c and C. Similarly,

$$
k v_h k_{0;[e]}^2 > k S_d(v_h) k_{0;j[e]} k (S_d(v_h) v_h) k_{0;j[e]}^2 > \frac{3}{4} k S_d(v_h) k_{0;j[e]}^2 \quad \text{Ch} j v_h j_{1;T_e}^2 : \quad (48)
$$

Putting (47) and (48) into (46) and summing over e ~ h;D we arrive at

$$
\begin{equation}
\begin{aligned}
& X \quad h S_d (r g_h) \mathbf{n} ; S_d (v_h) j_e i_e > \frac{e_1}{h} k S_d (v_h) k_{0; \sim h; D}^2 e_2 j v_h j_{1; \sim h}^2 \\
& e \stackrel{\sim}{=} \mathbf{e}_{h;D}
\end{aligned}
\tag{49}
\end{equation}
$$

with $e_1 = c_D c_1 = 2$ and $e_2$ regrouping the constants in (46), (47) and (48).

We shall also need to bound some norms of $g_h$ by those of $v_h$. To this end, we start from the following bound, analogous to equation (14b) in the proof of [28, Lemma 1]:

$$
jg_h j_1^2 + \frac{1}{h} \mathbf{e}^{\sim}_{h;D} \mathbf{X} \quad k S_d (g_h) k_{0;[e]}^2 6 \frac{c_3}{h} \mathbf{e}^{\sim}_{h;D} \mathbf{X} \quad k v_h k_{0;[e]}^2 : \quad (50)
$$

The original version of this bound is slightly different and does not include the factors $j_e$, nor the shift operator $S_d$. The present version is easily proven by scaling, having in mind Assumption 2. Note that

$$
\frac{1}{h} k \llbracket g_h \rrbracket k_{0;E \text{ ext}; \circ}^2 6 c_4 j g_h j_1^2 : \quad (51)
$$

Indeed, on any exterior (ctitious) edge $e_a$ adjacent to a boundary node a, so that $e_a$ is shared by extensions of interior mesh triangles, say $T_1; T_2$, we have

$$
jr g_h j_{T_1} \quad r \quad g_h j_{T_2} \quad j \quad 6 \quad \frac{c_4}{h} j g_h j_{1;T_1}^2 [T_2]
$$

with an h-independent constant $c_4$. Since $\llbracket g_h \rrbracket (a) = 0$, we can estimate

$$
k \llbracket g_h \rrbracket k_{e_a}^2 6 h^2 k r g_h j_{T_1} \quad r \quad g_h j_{T_2} k_{e_a}^2 6 c_4 h j g_h j_{1;T_1}^2 [T_2] :
$$
Introducing, on any edge $e$, the averages $f$ and $v_h$ of $f$ and $v_h$ on that edge, we have that
$e v_h = v_h jej$, and we continue the above calculations as

$$
\begin{align}
a_h(v_h; u \quad u_h) &= \underbrace{X}_{e \sim_h} (v_h \quad v_h; f)_{T_e^{\text{ext}}} + (v_h; f \quad f)_{T_e^{\text{ext}}} + (v_h; f)_{T_e^{\text{ext}}} \quad h \quad v_h \quad v_h; f \quad i_e \quad H_e \quad v_h f \, j T_e^{\text{ext}} \, j \nonumber \\
&= \underbrace{X}_{e \sim_h} (v_h \quad v_h; f)_{T_e^{\text{ext}}} + (v_h; f \quad f)_{T_e^{\text{ext}}} \quad h \quad v_h \quad v_h; f \quad i_e \quad H_e \nonumber \\
&\qquad 6 \text{ Chjv}_h j_1; \text{ kf } k_0; \quad + \text{ Chkv}_h k_0; \quad \text{jf } j_1; \quad + \text{ Ch}^2 \text{jv}_h j_1; \sim_h \text{ kf } k_{\sim_h} : \tag{54}
\end{align}
$$

To derive the last bound, we have used the following inequalities involving any boundary edge $e \sim_h$, shared by the element $T_e \in \{2, \dots, 6\}$ and the (curved) quadrilateral $T_e^{\text{ext}}$ (see Fig. 3 for an illustration):

$$
kv_h \quad v_h k_{T_e^{\text{ext}}} \quad 6 \quad \text{Chjv}_h j_1;_{T_e} ; \quad kv_h k_{T_e^{\text{ext}}} \quad 6 \quad Ck_{v_h} k_{T_e} ; \qquad (55)
$$

$$
\text{kf} \quad f \quad k_{T_e^{\text{ext}}} \quad 6 \quad \text{Chj} f \, j_1;_{T_e} \{ T_e^{\text{ext}} \} ; \quad kv_h \quad v_h k_e \quad 6 \quad \text{Chj} v_h j_1;_e : \qquad (56)
$$

All of these can be easily proven by scaling arguments and Poincaré-type inequalities. We have also used the fact that $H_e$ 6 Ch. Now, adding to these the trace inverse inequality $jv_hj_1;_{\sim_h} 6 \frac{C}{h}jv_hj_1;_{\sim_h}$ and the trace inequality $kf k_{\sim_h} 6 kf k_1;_{\sim_h}$ yields

$$
a_h(v_h; u \quad u_h) 6 Chkv_h k_1; \quad kf k_1; \quad 6 Ch 9 v_h 9_h \quad kf k_1; : \qquad (57)
$$

The passage from the full $H^1$ norm of $v_h$ to its triple norm in the last line, is justified by the Poincaré-type inequality (valid since $D$ is assumed of positive measure): for any $v \in 2H^1(\Omega)$

$$
kv k_0; \quad C_p \quad kr v k_0; \quad + j \quad D j \quad 1=2 \quad kv k_0; \quad D :
$$

From this and the fact that the boundary term on $\sim_{h,D}$ in the triple norm (42) can be interpreted as the
norm on $D$ when applied to $v_h \in \nabla_h^{\text{ext}}(\Omega)$, we deduce that

$$
kv_h k_0; \quad C_p \quad kr v_h k_0; \quad + h \quad 1=2 \quad kv_h k_{0; \sim_{h,D}} \quad C_p \quad 9 v_h 9_h : \qquad (58)
$$

Introducing the interpolant $I_h u \in \mathcal{H}_0^1(\Omega)$ (constructed by nodal interpolation over the mesh $T_h$ inside $\sim_h$), we derive from (57) using either the coercivity of the form $a_h$ from Lemma 1, or the inf-sup property from Lemma 2, depending on the choice of parameters $\cdot, \cdot$:

$$
9 \, I_h u \quad u_h \, 9_h \quad \sup_{w_h \in \mathcal{W}_h} \frac{a_h(w_h; I_h u \quad u_h)}{9 w_h \, 9_h} \quad \sup_{w_h \in \mathcal{W}_h} \frac{a_h(w_h; u \quad I_h u)}{9 w_h \, 9_h} + \text{Chkf } k_1; \quad : \qquad (59)
$$

Using the continuity of the form $a_h$ in the triple norm (evident by Cauchy-Schwarz and scaling), we arrive
at:

$$
9 I_h u \quad u_h 9_h \quad 6 \quad \text{Chkf } k_1; \quad + \quad C \quad 9 \quad u \quad I_h u \quad 9_h :
$$

By the Bramble Hilbert lemma [11, 12] and scaling arguments, we have

$$
9 u \quad I_h u 9_h \quad 6 \quad C \quad \text{ju} \quad I_h u j_1; \quad 6 \quad \text{Chju} j_2; :
$$

This allows us to obtain the announced error estimate, thanks to the triangle inequality

$$
9u \quad u_h 9_h \quad 6 \quad 9u \quad I_h u 9_h + 9 I_h u \quad u_h 9_h \quad 6 \quad \text{Chkf } k_1; + C \quad 9u \quad I_h u 9_h \quad 6 \quad \text{Ch}(juj_2; + \text{kf } k_1; ) : \quad (60)
$$
We now turn to the $L^2$ error estimate, and we initially focus on the symmetric interior penalty method (i.e. for $h=1$ and suitably large). This variant possesses the adjoint consistency property required for optimal convergence. To treat the errors introduced by the approximations in the gap $n \sim_h h$, we enlarge it to the set $h^{\text{fat}}$, which contains $n \sim_h h$ together with all elements $T^2 2 T_h$ adjacent to $h^{\sim_h}$. Note that $h^{\text{fat}}$ is a thin layer of width of order $h$, a fact that will be used via the following lemmas:

**Lemma 3.** For any $u \in L^1(\Omega)$

$$ \kappa u_{k_0; \frac{h}{\text{fat}}} = 6 C^p \bar{h} \kappa u_{k_1; \frac{h}{\text{fat}}} : $$

This lemma is proved in [37, Lemma 4.10]. The idea of the proof is to cover the band $h^{\text{fat}}$ by curves $S = f \times 2 : \text{dist}(x; \partial\Omega) = g$ for $2(0; \cdot)$ with of order $h$, apply the trace inequality on each $S$, and then integrate on $\Omega$.

**Lemma 4.** For any $u_h \in L^2(\Omega)$

$$ \kappa u_h k_{0; \frac{h}{\text{fat}}} = 6 C^p \bar{h} (\kappa u_h k_{0; \frac{h}{\text{fat}}} + \text{jump}_h) $$

with

$$ \text{jump}_h = \sum_{T^2 2 T_h} \sum_{T_h^{\text{ext}}} \sum_{T_h^{\text{ext},0}} \sum_{T_h^{\text{ext},0}} \frac{1}{h} k[r u_h] k_{0;e}^2 : $$

**Proof**

The idea is to apply the preceding lemma to $r u_h$, but the difficulty is that $r u_h$ is not in $H^1(\Omega)$, because of the discontinuities that arise in the gap $n \sim_h h$. To circumvent this issue, we can construct a continuous counterpart $U_h$ of $r u_h$. For example, we can set $U_h = I_h r u_h$ where $I_h$ is a Cløment interpolation operator to $P^1$ functions. Then $kU_h k_0; 6 C \kappa u_h k_0;$, $kr u_h U_h k_0; 6 C \text{jump}_h$, and $kr U_h k_0; 6 C \text{jump}_h$; by scaling arguments. We conclude using Lemma 3:

$$ \kappa u_h k_{0; \frac{h}{\text{fat}}} = 6 C^p \bar{h} \kappa u_h U_h k_{0; \frac{h}{\text{fat}}} + \kappa U_h k_{0; \frac{h}{\text{fat}}} = 6 C^p \bar{h} U_h k_1; $$

With the interpolation estimates from above, this gives

$$ \kappa u_h k_{0; \frac{h}{\text{fat}}} = 6 C^p \bar{h} (\kappa u_h k_{0; \frac{h}{\text{fat}}} + \text{jump}_h) : $$

2

**Theorem 2** (optimal $L^2$-error estimate for the symmetric interior penalty method). Suppose $h=1$, suitably large, $f \in L^2(\Omega)$ and $u \in L^2(\Omega)$. The following $L^2$ error estimate holds then

$$ \text{jump}_h = 6 C^p \bar{h} (\text{jump}_h + \text{jump}_k) : $$

**Proof**

We proceed by the Aubin-Nitsche trick. Let $w \in L^2(\Omega)$ be the solution to

$$ w = u u_h; \text{ in } \Omega; \tag{61a} $$

$$ w = 0; \text{ on } \partial\Omega; \tag{61b} $$

$$ \text{jump}_h = 0; \text{ on } \partial\Omega; \tag{61c} $$

Then

$$ \text{jump}_h = a_h(u u_h; w) = a_h(w; u u_h); $$
by the symmetry of $a_h$ for the symmetric variant of the method with $= 1$. Take $w_h = I_h w$ $V_h$ as the nodal interpolant of $w$. Then, using the modified Galerkin orthogonality (53) we rewrite

$$ \text{ku } u_h k_{0;}^2 = a_h(w w_h; u u_h) + \sum_{e \sim_h} X (w_h; f)_{T_e^{\text{ext}}} h w_h; f i_e H_e : $$

Proceeding as in the proof of the $H^1$ error by introducing the averages $f$ and $w_h$ of $f$ and $w_h$ on every edge $e$, we continue the above calculations as

$$ \text{ku } u_h k_{0;}^2 = \sum_{6 \text{ 9 w } w_h} u_h w_h \sum_{9 \text{ u } u_h} u_h w_h + \sum_{e \sim_h} X (w_h w_h; f)_{T_e^{\text{ext}}} + (w_h; f f)_{T_e^{\text{ext}}} \\ \qquad \times (hw_h w_h; f i_e H_e) \\ \qquad \sum_{e \sim_h} 6 \text{ Ch}^2 \text{jw}_{j2;} (\text{juj}_{j2;} + \text{kf} k_1;) + \text{Ch}(\text{jw}_{h1;} \text{ fat} \text{ kf} k_0; \text{ fat} \text{ + kw}_h k_0; \text{ fat} \text{ jf} j_1; \text{ fat}) \\ \qquad + \text{Ch}^2 \text{jw}_{j1; \sim_h} \text{ kf} k_0; \text{ fat} : $$

Here, we have used the already proven error estimate (60), and bounds (55). To gain another $h$ in the terms multiplied by the first power of $h$, we use Lemmas 3 and 4 :

$$ \text{ku } u_h k_{0;}^2 = \sum_{6 \text{ Ch}^2 \text{jw}_{j2;}} (\text{juj}_{j2;} + \text{kf} k_1;) + \sum_{\text{Ch}^2((\text{jw}_{h1;} + \text{jw}_{h2;h;})\text{kf} k_1; + \text{jw}_{h1;} \text{kr} f k_1;)} \\ + \sum_{\text{Ch}^2 \text{jw}_{h1; \sim_h}} \text{kf} k_0; \text{ fat} : $$

By interpolation estimates, $\text{jw}_{h1; \sim_h} + \text{jw}_{h1;} + \text{jw}_{h2;h;}$ $6 \text{ Ckw}_2;$. By the trace inequality, $\text{kf} k_0; \sim_h$ $6 \text{ Ckf} k_1;$. Thus, the estimate above leads to

$$ \text{ku } u_h k_{0;}^2 = \sum_{6 \text{ Ch}^2(\text{juj}_2; + \text{kf} k_2;)} \text{kw}_2; $$

and we conclude by recalling $\text{kw}_2; 6 \text{ Cku } u_h k_0;$. $2$

Remark 8. For other choices of parameters, i.e. $1$ with sufficiently large, or $= 1; = 0$, we can prove the sub-optimal $L^2$ error estimate

$$ \text{ku } u_h k_0; = \sum_{6 \text{ Ch}^{3=2}(\text{juj}_2; + \text{kf} k_2;)} \text{Ch}^2; $$

This result should not be surprising, since it also holds for the unsymmetric Nitsche's method on body-tted grids. To this end, we proceed again by the Aubin-Nitsche trick, starting by introducing $w$ $2 H^2(\cdot)$ as in the proof of Theorem 2. This time, however, the bilinear form is not symmetric, and some adjustments need to be made. Then

$$ \text{ku } u_h k_{0;}^2 = a_h(u u_h; w) \\ a_h(w; u u_h) \quad (1) \quad h[\llbracket u u_h \rrbracket]; \text{ fr } w \mathfrak{g} i_{E^{\text{ext}; o}} \\ (1) \quad h S_d(r w) n; S_d(u u_h) j_e i_e \\ \sum_{e \sim_h; D} X (w_h; f)_{T_e^{\text{ext}}} h w_h; f i_e H_e + c; \quad (62) \\ = a_h(w w_h; u u_h) + \sum_{e \sim_h} (w_h; f)_{T_e^{\text{ext}}} h w_h; f i_e H_e + c; $$

where

$$ c = (1) h[\llbracket u u_h \rrbracket]; \text{ fr } w \mathfrak{g} i_{E^{\text{ext}; o}} \quad (1) \quad h S_d(r w) n; S_d(u u_h) j_e i_e; \\ \sum_{e \sim_h; D} X $$

Other than $c$, all the terms in the right-hand side of (62), can be treated as in the proof of Theorem 2.
For c, we proceed as follows

$$c 6 C \frac{p}{\bar{h}} 9 u \ u_h 9_h \ kr w n k_{0;\sim_{h,D}} + kr w n k_{0;E^{\text{ext}};o} :$$

To treat the contributions on the edges in $E^{\text{ext}};o$, we apply the following inverse inequality

$$kr w n k_{0;E^{\text{ext}};o} 6 C \frac{1}{\bar{h}} kr w k_{0;\bar{h};\text{fat}} + \frac{p}{\bar{h}} jr w j_{1;\bar{h};\text{fat}} 6 C kr w k_{1;\bar{h}} ;$$

derived applying Lemma 3, and thinking about the edges in $E^{\text{ext}};o$ as part of quadrilaterals in $T_h^{\sim \text{ext}}$. Collecting all these contributions and using the estimate for $9 u \ u_h 9_h$ gives

$$\text{ku } u_h k_{0;2}^2 + 6 \text{Ch}^2(\text{juj}_2; + \text{kf } k_2; )\text{kwk}_2; + \text{Ch}^{3=2}(\text{juj}_2; + \text{kf } k_1; ) kr w n k_{0;\sim_{h,D}} + kr w k_1;$$

and we conclude by recalling $\text{kwk}_2; 6 \text{Cku } u_h k_0; .$

### 5. Neumann Boundary Conditions for Linear Elasticity

In the numerical tests that follow, we also consider the equations of (compressible) isotropic linear elasticity. Their strong form is given as

$$r \quad (u) = b \quad \text{in } ; \qquad (63a)$$

$$u = u_D \quad \text{on } D;h; \qquad (63b)$$

$$n = t_N \quad \text{on } N;h; \qquad (63c)$$

where $u$ is the displacement field, $u_D$ its value on the Dirichlet boundary $D;h$, $t_N$ the normal traction along the Neumann boundary $N;h$, and $b$ a body force. We of course assume that $@ = \frac{D;h}{N;h}$ and $D;h \setminus N;h = ;$. The stress is a linear function of $u$, according to the constitutive model

$$(u) = 2 \quad (u) + (r \ u)I :$$

The proposed SBM variational form of (63) can be derived in a similar way to the case of the Poisson equation. Using the identities:

$$[! h] = !_h^+ n^+ + !_h n ; \qquad (64)$$

$$[v_h] = v_h^+ n^+ + v_h n ; \qquad (65)$$

$$[!_h v_h] = f!_h g : [[v_h]] + [[!_h]] f v_h g; \qquad (66)$$

where $!_h$ is a second-order tensor function and $v_h$ is a vector function. Hence we have:

$$ \begin{aligned} & ("(w_h); 2"(u_h))_{\sim_h} + ("(w_h); 2"(u_h))_{T_h^{\sim \text{ext}}} + (r w_h; r u_h)_{\sim_h} + (r w_h; r u_h)_{T_h^{\sim \text{ext}}} \\ & h [[w_h]] ; f r u_h I + 2"(u_h)g i_{E^{\text{ext}};o} \quad h f r w_h I + 2"(w_h)g ; [[u_h]] i_{E^{\text{ext}};o} \\ & h w_h ; (r u_h I + 2"(u_h)) n i_{D;h} \quad h(r w_h I + 2"(w_h)) n ; u u_D i_{D;h} \\ & \qquad + h h {}^1[[w_h]] ; [[u_h]] i_{E^{\text{ext}};o} + h h {}^1 w_h ; u_h u_D i_{D;h} \\ & = (w_h ; b)_{\sim_h} + (w_h ; b)_{T_h^{\sim \text{ext}}} + h w_h ; t_N i_{N;h} : \end{aligned} \qquad (67) $$

As in the case of the Poisson equation, the variational form 67 can be approximated as

$$ ("(w_h); 2"(u_h))_{\sim_h} + (r w_h; r u_h)_{\sim_h} + \begin{matrix} X \\ e \end{matrix} h"(w_h); 2"(u_h)H_e i_e + h r w_h; r u_h H_e i_e $$
$$
\begin{align*}
& X \frac{jd_{\alpha}j}{2} [[r w_h d_{\alpha}]_{;\alpha} : f r u_h I + 2 \quad (u_h)g_{\alpha} + f r w_h I + 2 \quad (w_h)g_{\alpha} : [[r u_h d_{\alpha}]_{;\alpha}} \\
& \quad \alpha 2 N (\sim_h) \\
& X hS_d(w_h) ; ((r u_h I + 2 \quad (w_h))n)j_e i_e + h(r u_h I + 2 \quad (w_h))n ; (S_d(u_h) u_D)j_e i_e \\
& \quad e^{\sim_h D} \\
& X + \frac{4h}{\alpha 2 N (\sim_h)} jd_{\alpha}j [[r w_h d_{\alpha}]_{;\alpha} [[r u_h d_{\alpha}]_{;\alpha}} + X \quad e^{\sim_h D} \quad X \\
& \quad hS_d(w_h) ; (S_d(u_h) u_D)j_e i_e \\
& = (w_h ; b)_{\sim_h} + X \quad e^{\sim_h} \quad h w_h ; b H_e i_e + X \quad e^{\sim_h N} \quad hS_d(w_h) ; h_N(M_h(\star))j_e i_e : \quad (68)
\end{align*}
$$

<figure><img src="image_5.png" /><figcaption>Figure 4: General problem schematic for the active domain (grey), true boundary (blue), surrogate boundary (red)</figcaption></figure>

6. Numerical results

We present the results from a series of two-dimensional numerical experiments that demonstrate the
theoretical ndings of Section 4. Our general approach relies on the method of manufactured solutions and
involves embedding geometries (both analytic and polygonal) on a series of meshes of increased re nement.
A schematic representation of the experimental setup is provided in Figure 4, depicting an embedded circular
shape on a unit square domain.

More speci cally, since embedded geometries may arbitrarily intersect the grids, we also apply speci c
rotations to the latter (with respect to the embedded shapes) to test the e ect on the numerical results of
such perturbations. In fact, the newly conceptualized method requires computation of additional geometric
quantities on the surrogate boundary (with respect to a traditional SBM). Grid rotations provide a means
to examine a high number of extension arrangements, i.e. robustness. The convergence of the L²-norm and
the H¹-seminorm of the error are assessed, along with the condition number, for both the Poisson and linear
elasticity equations.

The primary motivator for the proposed method was, by and large, the development of Neumann bound-
ary conditions that guarantee optimal convergence rates. Although a majority of the experiments involve
only embedded Neumann boundaries, the array of tests was expanded to also include Dirichlet boundary
conditions. The numerical experiments on the Poisson equation encompass both the symmetric and anti-
symmetric Nitsche formulations for the weak enforcement of Dirichlet conditions, along with the inclusion of
a patch test. Likewise, a simple bending beam test was performed for linear elasticity, which included both
homogeneous displacement and homogeneous traction boundary conditions. Furthermore, optimal conver-
gence was also achieved with quadrilateral elements, demonstrating the exibility of the method beyond
standard triangular nite elements.

6.1. Patch test for the Poisson problem

A patch test experiment was performed to assess the ability of the method to match an a ne exact
solution. In general, passing a patch test does not guarantee convergence nor stability of a numerical method.
Yet, it is an important sanity check for the proposed conceptualization of the SBM, which involves solution
extensions and approximate integration over the gap region. A ne solutions are relevant in engineering
applications, since they imply a constant ux scenario in the case of the Poisson problem or a constant
strain scenario in the case of linear elasticity.

The proposed SBM variant possesses the partition of unity property and passes the patch test, since
a ne solutions can be exactly represented in the gap region. We considered the Poisson problem with
three simple geometries (circle, rotated square, and star) embedded on a unit square domain with a regular
background mesh, as seen in Figure 5. The solutions are u = x + y, u = x, and u = y for the circle, rotated
square, and star respectively. A Dirichlet boundary condition is strongly enforced on the outer perimeter
of the unit square, while the shifted Neumann boundary condition is applied to the surrogate boundary in
red. A visualization of the surrogate extensions are provided in Figure 6, along with contours of the nodal
error between the approximate and exact solutions u_h u_exact. It is easily seen that the numerical error is
within machine precision.

<figure><img src="image_6.png" /><figcaption>Figure 5: Patch test solutions u_h on surrogate domains ~_h. (blue) and ~_h (red)</figcaption></figure>

<figure><img src="image_7.png" /><figcaption>Figure 6: Patch test results with u_h u_exact contours. (blue) and T<sup>ext</sup> extensions (red)</figcaption></figure>
## 6.2. Poisson problem with manufactured solutions

The aim of these numerical experiments is to analyze convergence properties of the method for a high number of extension arrangements. Rotating the background mesh, while xing the embedded geometry and manufactured solution, isolates e ects imparted on the quality of the solution by the geometry of the element extensions $T_e^{ext}$. Both a smooth, analytical shape (circle) and a concave, polygonal shape (star) were immersed into triangular background meshes. Nine increments of grid rotation (from zero to 45 degrees) and seven levels of grid re nement were applied. Results from a boundary- tted, primal formulation are also included for comparison.

The geometric setup of the first test involves a circular boundary of radius 0.25 and centered at [0.6,0.5]. The computational grids are unstructured triangular meshes rotated around [0.5,0.5]: this o set ensures variability in the cuts for each rotation. Neumann boundary conditions are applied along the surrogate boundary associated with the circular shape. The second test involves a star-shaped polygon (ve-point star) centered at [0.5,0.5] and immersed into a structured background triangular mesh. Neumann boundary conditions are applied on the surrogate boundary along the star shape for $x > 0.5$ and Dirichlet conditions for $x \le 0.5$. As in the patch test, Dirichlet conditions are enforced strongly on the outer perimeter for both tests. The analytical solution and corresponding forcing function are

$$u(x; y) = \sin(4x)\sin(4y); \qquad (69)$$

$$f(x; y) = 32^2 \sin(4x)\sin(4y); \qquad (70)$$

which were deduced from the strong form of the Poisson equation, using the method of manufactured solutions. Boundary conditions are specified accordingly.

For visualization purposes, a sampling of computed solutions from various rotations are included in Figures 7 and 9. The results displayed in Figures 8 and 10 show that the convergence rates of the L²-norm and the H¹-seminorm are optimal. In terms of the condition number (A) associated with the algebraic problem, we see that the proposed method maintains the expected scaling of (A) h². In Figure 10, both the symmetric ($= 1$, $= 10$) and anti-symmetric ($= 1$, $= 0$) Nitsche formulations are simulated and compared to the primal, boundary- tted case. The symmetric Nitsche formulation seems more accurate in the L²-norm of the error, but at the expense of higher condition numbers.

<figure><img src="image_8.png" /><figcaption>Figure 7: Visualizations of $u^h$ for various degrees of rotation, superimposed on background grid.</figcaption></figure>

Similar convergence tests with the same manufactured solution were performed using quadrilateral nite elements. For boundaries, a square with a side length of $1 = 0:48$ and centered at [0:5; 0:5] and a concave ower-like geometry are considered. The coordinates of the ower-like boundary are parametrized as functions of the angle $\langle : $

$$ \begin{aligned} (x(\langle)) &= 0:5 + (0:05 + 0:24 \sin(7 \langle)) \cos(\langle)); \\ y(\langle) &= 0:5 + (0:05 + 0:24 \sin(7 \langle)) \sin(\langle); \end{aligned} \qquad (71) $$
<figure><img src="image_9.png" /><figcaption>(a) L²-norm for all rotations.</figcaption></figure>

<figure><img src="image_10.png" /><figcaption>(b) H¹-seminorm for all rotations.</figcaption></figure>

<figure><img src="image_11.png" /></figure>
obtained with triangular grids.

<figure><img src="image_12.png" /><figcaption>(a) and ~ for rotation of 40° (b) $u^h$ for rotation of 40° (c) and ~ for rotation of 0° (d) $u^h$ for rotation of 0°

Figure 11: True boundary (blue), surrogate boundary (~-red) and $u^h$ on the surrogate domain ~$^h$.</figcaption></figure>

<figure><img src="image_13.png" /><figcaption>Figure 12: Convergence rates and condition numbers for the Poisson problem with a square boundary over Cartesian grids. Because the results for body-tted grids are typically clustered more tightly, we do not perform grid rotations for this case.</figcaption></figure>

<figure><img src="image_14.png" /></figure>
boundary. As before, the computational grids were incrementally rotated from 0 degrees to 45 degrees for seven levels of grid re nement. The analytical solution was chosen to be

$$
\begin{align}
& ( u_x(x; y) = \sin(2x) \sin(2y); \notag \\
& u_y(x; y) = \cos(2x) \cos(2y); \tag{72}
\end{align}
$$

<figure><img src="image_15.png" /><figcaption>(a) L²-norm for all rotations.</figcaption></figure>

<figure><img src="image_16.png" /><figcaption>(b) H¹-seminorm for all rotations.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>L²-norm and H¹-seminorm</figcaption></figure>
<figure><img src="image_18.png" /><figcaption>Figure 16: Convergence rates for the linear elasticity problem with square boundary and Cartesian grids.</figcaption></figure>

<figure><img src="image_19.png" /></figure>
<figure><img src="image_20.png" /><figcaption>Figure 18: Cantilever beam bending test: Surrogate domain ~h with true (blue) and surrogate boundaries ~h (green, red), and deformed con guration.</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>Figure 19: Cantilever beam bending test: tip vertical de ection for the Boundary tted method (red) and SBM (blue). The reference solution of the Euler Beam Theory is in black. n is the number of boundary segments.</figcaption></figure>

of the integration procedure on the gap relies on the concept of a distance, and because extensions of the solution in the form of Taylor expansions are used to evaluate the solution in the gap. Hence, no cut-cell integration procedure is performed. A series of numerical experiments proved the consistency, stability, robustness, and optimal accuracy of the proposed approach.

Acknowledgments

G. Scovazzi has been partially supported by the National Science Foundation (Division of Mathematical Sciences), with Grant DMS 2207164 and Grant DMS 2409919.
References

[1] Nabil Atallah, Claudio Canuto, and Guglielmo Scovazzi. Analysis of the Shifted Boundary Method for the Poisson problem in domains with corners. Mathematics of Computation, 90(331):2041 2069, 2021.

[2] Nabil M Atallah, Claudio Canuto, and Guglielmo Scovazzi. Analysis of the Shifted Boundary Method for the Stokes problem. Computer Methods in Applied Mechanics and Engineering, 358:112609, 2020.

[3] Nabil M Atallah, Claudio Canuto, and Guglielmo Scovazzi. The Shifted Boundary Method for solid mechanics. International Journal for Numerical Methods in Engineering, 122(20):5935 5970, 2021.

[4] Nabil M. Atallah, Claudio Canuto, and Guglielmo Scovazzi. The shifted boundary method for solid mechanics. International Journal for Numerical Methods in Engineering, 122(20):5935 5970, 2021.

[5] Nabil M Atallah, Claudio Canuto, and Guglielmo Scovazzi. The high-order Shifted Boundary Method and its analysis. Computer Methods in Applied Mechanics and Engineering, 394:114885, 2022.

[6] Santiago Badia, Francesc Verdugo, and Alberto F Martín. The aggregated untted nite element method for elliptic problems. Computer Methods in Applied Mechanics and Engineering, 336:533 553, 2018.

[7] Silvia Bertoluzza, Mourad Ismail, and Bertrand Maury. The Fat Boundary Method: Semi-discrete scheme and some numerical experiments. In Domain decomposition methods in science and engineering, pages 513 520. Springer, 2005.

[8] Silvia Bertoluzza, Mourad Ismail, and Bertrand Maury. Analysis of the fully discrete Fat Boundary Method. Numerische Mathematik, 118(1):49 77, 2011.

[9] Daniele Bo and Lucia Gastaldi. A nite element approach for the Immersed Boundary Method. Computers & structures, 81(8):491 501, 2003.

[10] James H Bramble, Todd Dupont, and Vidar Thomøe. Projection methods for Dirichlet's problem in approximating polygonal domains with boundary-value corrections. Mathematics of Computation, 26(120):869 879, 1972.

[11] James H Bramble and SR Hilbert. Estimation of linear functionals on sobolev spaces with application to fourier transforms and spline interpolation. SIAM Journal on Numerical Analysis, 7(1):112 124, 1970.

[12] James H Bramble and SR Hilbert. Bounds for a class of linear functionals with applications to Hermite interpolation. Numerische Mathematik, 16(4):362 369, 1971.

[13] James H Bramble and J Thomas King. A robust nite element method for nonhomogeneous Dirichlet problems in domains with curved boundaries. mathematics of computation, 63(207):1 17, 1994.

[14] James H Bramble and J Thomas King. A nite element method for interface problems in domains with smooth boundaries and interfaces. Advances in Computational Mathematics, 6(1):109 138, 1996.

[15] Erik Burman. Ghost penalty. Comptes Rendus Mathematique, 348(21-22):1217 1220, 2010.

[16] Erik Burman, Susanne Claus, Peter Hansbo, Mats G Larson, and AndrØ Massing. CutFEM: Discretizing geometry and partial di erential equations. International Journal for Numerical Methods in Engineering, 104(7):472 501, 2015.

[17] Erik Burman, Daniel Elfverson, Peter Hansbo, Mats G Larson, and Karl Larsson. Shape optimization using the Cut Finite Element Method. Computer Methods in Applied Mechanics and Engineering, 328:242 261, 2018.

[18] Erik Burman and Miguel A FernÆndez. An untted Nitsche method for incompressible uid structure interaction using overlapping meshes. Computer Methods in Applied Mechanics and Engineering, 279:497 514, 2014.

[19] Erik Burman and Peter Hansbo. Fictitious domain nite element methods using cut elements: I. A stabilized Lagrange multiplier method. Computer Methods in Applied Mechanics and Engineering, 199(41-44):2680 2686, 2010.

[20] Erik Burman and Peter Hansbo. Fictitious domain nite element methods using cut elements: II. A stabilized Nitsche method. Applied Numerical Mathematics, 62(4):328 341, 2012.

[21] Erik Burman and Peter Hansbo. Fictitious domain methods using cut elements: Iii. a stabilized Nitsche method for Stokes' problem. ESAIM: Mathematical Modelling and Numerical Analysis, 48(3):859 874, 2014.

[22] Erik Burman, Peter Hansbo, and Mats Larson. A cut nite element method with boundary value correction. Mathematics of Computation, 87(310):633 657, 2018.

[23] Erik Burman, Peter Hansbo, and Mats G Larson. A cut nite element method with boundary value correction for the incompressible Stokes equations. In European Conference on Numerical Mathematics and Advanced Applications, pages 183 192. Springer, 2017.

[24] Erik Burman, Peter Hansbo, and Mats G Larson. Dirichlet boundary value correction using Lagrange multipliers. arXiv preprint arXiv:1903.07104, 2019.

[25] Bernardo Cockburn, Weifeng Qiu, and Manuel Solano. A priori error analysis for HDG methods using extensions from subdomains to achieve boundary conformity. Mathematics of Computation, 83(286):665 699, 2014.

[26] Bernardo Cockburn and Manuel Solano. Solving Dirichlet boundary-value problems on curved domains by extensions from subdomains. SIAM Journal on Scientific Computing, 34(1):A497 A519, 2012.

[27] Bernardo Cockburn and Manuel Solano. Solving convection-di usion problems on curved domains by extensions from subdomains. Journal of Scientific Computing, 59(2):512 543, 2014.

[28] J Haydel Collins, Alexei Lozinski, and Guglielmo Scovazzi. A penalty-free shifted boundary method of arbitrary order. Computer Methods in Applied Mechanics and Engineering, 417:116301, 2023.

[29] Oriol ColomØs, Alex Main, LØo Nouveau, and Guglielmo Scovazzi. A weighted Shifted Boundary Method for free surface ow problems. Journal of Computational Physics, 424:109837, 2021.

[30] StØphane Cotin, Michel Duprez, Vanessa Lleras, Alexei Lozinski, and Killian Vuillemot. -FEM: an e cient simulation tool using simple meshes for problems in structure mechanics and heat transfer, 2022.

[31] Michel Duprez, Vanessa Lleras, and Alexei Lozinski. A new -FEM approach for problems with natural boundary conditions. Numerical Methods for Partial Di erential Equations, 39(1):281 303, 2023.

[32] Michel Duprez, Vanessa Lleras, and Alexei Lozinski. -FEM: an optimally convergent and easily implementable immersed
boundary method for particulate ows and stokes equations.
57(3):1111 1142, 2023.

[33] Michel Duprez, Vanessa Lleras, Alexei Lozinski, and Killian Vuillemot. An immersed boundary method by approach to solve the heat equation, 2022. -FEM

[34] Michel Duprez, Vanessa Lleras, Alexei Lozinski, and Killian Vuillemot. phi-FEM for the heat equation: optimal convergence on untted meshes in space. arXiv preprint arXiv:2303.12013, 2023.

[35] Michel Duprez and Alexei Lozinski. -FEM: a nite element method on domains de ned by level-sets. SIAM Journal on Numerical Analysis, 58(2):1008 1028, 2020.

[36] Alexander D ster, Jamshid Parvizian, Zhengxiong Yang, and Ernst Rank. The Finite Cell Method for three-dimensional problems of solid mechanics. Computer methods in applied mechanics and engineering, 197(45):3768 3782, 2008.

[37] Charles M Elliott and Thomas Ranner. Finite element analysis for a coupled bulk surface partial di erential equation. IMA Journal of Numerical Analysis, 33(2):377 402, 2013.

[38] Roland Glowinski, Tsorng-Whay Pan, and Jacques Periaux. A citious domain method for Dirichlet problem and applications. Computer Methods in Applied Mechanics and Engineering, 111(3-4):283 303, 1994.

[39] Anita Hansbo and Peter Hansbo. An untted nite element method, based on Nitsche's method, for elliptic interface problems. Computer methods in applied mechanics and engineering, 191(47):5537 5552, 2002.

[40] Klaus Hllig. Finite element methods with B-splines. SIAM, Philadelphia, 2003.

[41] Klaus Hllig, Ulrich Reif, and Joachim Wipper. Weighted extended B-spline approximation of Dirichlet problems. SIAM Journal on Numerical Analysis, 39(2):442 462, 2001.

[42] David Kamensky, Ming-Chen Hsu, Yue Yu, John A Evans, Michael S Sacks, and Thomas JR Hughes. Immersogeometric cardiovascular uid structure interaction analysis with divergence-conforming B-splines. Computer Methods in Applied Mechanics and Engineering, 314:408 472, 2017.

[43] Kangan Li, Nabil M Atallah, G Alex Main, and Guglielmo Scovazzi. The Shifted Interface Method: a exible approach to embedded interface computations. International Journal for Numerical Methods in Engineering, 121(3):492 518, 2020.

[44] Kangan Li, Nabil M Atallah, Antonio Rodr guez-Ferran, Dakshina M Valiveti, and Guglielmo Scovazzi. The Shifted Fracture Method. International Journal for Numerical Methods in Engineering, 122(22):6641 6679, 2021.

[45] Kangan Li, Andrea Gorgi, Riccardo Rossi, and Guglielmo Scovazzi. The shifted boundary method for contact problems. Computer Methods in Applied Mechanics and Engineering, 440:117940, 2025.

[46] Kangan Li, Antonio Rodr guez-Ferran, and Guglielmo Scovazzi. A blended Shifted-Fracture/Phase-Field framework for sharp/di use crack modeling. International Journal for Numerical Methods in Engineering, 124(4):998 1030, 2023.

[47] Kangan Li, Antonio Rodr guez-Ferran, and Guglielmo Scovazzi. The simple Shifted Fracture Method. International Journal for Numerical Methods in Engineering, 124:2837 2875, 2023.

[48] Chuanqi Liu and WaiChing Sun. Shift boundary material point method: an image-to-simulation work ow for solids of complex geometries undergoing large deformation. Computational Particle Mechanics, 7:291 308, 2020.

[49] Alexei Lozinski. A new citious domain method: Optimal convergence without cut elements. Comptes Rendus Mathematique, 354(7):741 746, 2016.

[50] Alexei Lozinski. CutFEM without cutting the mesh cells: A new way to impose Dirichlet and Neumann boundary conditions on untted meshes. Computer Methods in Applied Mechanics and Engineering, 356:75 100, 2019.

[51] Alex Main and Guglielmo Scovazzi. The Shifted Boundary Method for embedded domain computations. Part I: Poisson and Stokes problems. Journal of Computational Physics, 372:972 995, 2018.

[52] Alex Main and Guglielmo Scovazzi. The Shifted Boundary Method for embedded domain computations. Part II: Linear advection di usion and incompressible Navier Stokes equations. Journal of Computational Physics, 372:996 1026, 2018.

[53] AndrØ Massing, Mats Larson, Anders Logg, and Marie Rognes. A Nitsche-based cut nite element method for a uid-structure interaction problem. Communications in Applied Mathematics and Computational Science, 10(2):97 120, 2015.

[54] Jamshid Parvizian, Alexander D ster, and Ernst Rank. Finite Cell Method. Computational Mechanics, 41(1):121 133, 2007.

[55] T R berg and F Cirak. Subdivision-stabilised immersed B-spline nite elements for moving boundary ows. Computer Methods in Applied Mechanics and Engineering, 209:266 283, 2012.

[56] T R berg and F Cirak. A xed-grid B-spline nite element technique for uid structure interaction. International Journal for Numerical Methods in Fluids, 74(9):623 660, 2014.

[57] B Schott, U Rasthofer, V Gravemeier, and WA Wall. A face-oriented stabilized Nitsche-type extended variational multiscale method for incompressible two-phase ow. International Journal for Numerical Methods in Engineering, 104(7):721 748, 2015.

[58] Ting Song, Alex Main, Guglielmo Scovazzi, and Mario Ricchiuto. The Shifted Boundary Method for hyperbolic systems: Embedded domain computations of linear waves and shallow water ows. Journal of Computational Physics, 369:45 79, 2018.

[59] Jens Visbech, Allan P Engsig-Karup, and Mario Ricchiuto. A spectral element solution of the poisson equation with shifted boundary polynomial corrections: in uence of the surrogate to true boundary mapping and an asymptotically preserving robin formulation. Journal of Scientific Computing, 102(1):11, 2025.

[60] Fei Xu, Dominik Schillinger, David Kamensky, Vasco Varduhn, Chenglong Wang, and Ming-Chen Hsu. The tetrahedral Finite Cell Method for uids: Immersogeometric analysis of turbulent ow around complex geometries. Computers & Fluids, 141:135 154, 2016.

[61] Lucy Zhang, Axel Gerstenberger, Xiaodong Wang, and Wing Kam Liu. Immersed nite element method. Methods in Applied Mechanics and Engineering, 193(21-22):2051 2067, 2004.