# A spectral element solution of the Poisson equation with shifted boundary polynomial corrections: influence of the surrogate to true boundary mapping and an asymptotically preserving Robin formulation

Jens Visbech, Allan Engsig-Karup, Mario Ricchiuto

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

## To cite this version:

Jens Visbech, Allan Engsig-Karup, Mario Ricchiuto. A spectral element solution of the Poisson equation with shifted boundary polynomial corrections: influence of the surrogate to true boundary mapping and an asymptotically preserving Robin formulation. Journal of Scientific Computing, 2025, 102 (1), pp.11. 10.1007/s10915-024-02713-z. hal-04829220

HAL Id: hal-04829220

https://inria.hal.science/hal-04829220v1

Submitted on 10 Dec 2024

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published 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_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
# A SPECTRAL ELEMENT SOLUTION OF THE POISSON EQUATION WITH SHIFTED BOUNDARY POLYNOMIAL CORRECTIONS: INFLUENCE OF THE SURROGATE TO TRUE BOUNDARY MAPPING AND AN ASYMPTOTICALLY PRESERVING ROBIN FORMULATION

AN ARXIV PREPRINT

● **Jens Visbech**

Department of Applied Mathematics and Computer Science
Technical University of Denmark
Kongens Lyngby, 2800, Denmark
jvis@dtu.dk

● **Allan Peter Engsig-Karup**

Department of Applied Mathematics and Computer Science
Technical University of Denmark
Kongens Lyngby, 2800, Denmark
apek@dtu.dk

● **Mario Ricchiuto**

INRIA, Univ. Bordeaux, CNRS, Bordeaux INP, IMB, UMR 5251,
200 Avenue de la Vieille Tour, 33405 Talence cedex, France
mario.ricchiuto@inria.fr

October 27, 2023

## ABSTRACT

We present a new high-order accurate spectral element solution to the two-dimensional scalar Poisson equation subject to a general Robin boundary condition. The solution is based on a simplified version of the shifted boundary method employing a continuous arbitrary order *hp*-Galerkin spectral element method as the numerical discretization procedure. The simplification relies on a polynomial correction to avoid explicitly evaluating high-order partial derivatives from the Taylor series expansion, which traditionally have been used within the shifted boundary method. In this setting, we apply an extrapolation and novel interpolation approach to project the basis functions from the *true* domain onto the approximate *surrogate* domain. The resulting solution provides a method that naturally incorporates curved geometrical features of the domain, overcomes complex and cumbersome mesh generation, and avoids problems with small-cut-cells. Dirichlet, Neumann, and general Robin boundary conditions are enforced weakly through: i) a generalized Nitsche's method and ii) a generalized Aubin's method. For this, a consistent asymptotic preserving formulation of the embedded Robin formulations is presented.

We present several numerical experiments and analysis of the algorithmic properties of the different weak formulations. With this, we include convergence studies under polynomial, *p*, increase of the basis functions, mesh, *h*, refinement, and matrix conditioning to highlight the spectral and algebraic convergence features, respectively. This is done to assess the influence of errors across variational formulations, polynomial order, mesh size, and mappings between the true and surrogate boundaries.
**Keywords** Spectral element method · shifted boundary method · high-order numerical method · embedded methods · Poisson problem · elliptic problem · Dirichlet, Neumann, and Robin boundary conditions.

# 1 Introduction

The need for accurate representation of solutions to partial differential equations (PDEs) – that take complex geometry into account – appears frequently in real-world science and engineering applications. Some applications largely involve smooth fields, for example, free surface dispersive wave propagation, electro-magnetics, and electrostatics, among others. In these cases, very precise predictions can be obtained using, e.g., high-order finite elements on conformal meshes. However, the computational time spent meshing is often relatively high. The presence of moving boundaries, e.g., as in free surface potential flow wave models, further complicates this challenging task. Unfitted techniques – such as immersed and embedded methods – allow to work with meshes agnostic of the domain boundaries, hereby lighting the meshing task and offering greater geometrical flexibility. Another important aspect is the type of boundary condition being imposed. Many applications, such as some of those mentioned above, involve a mix of Dirichlet, Neumann, as well as Robin conditions. The latter plays an important role in many problems. To mention a few: relatively simple elliptic and parabolic PDEs, as heat transfer and fluid flow [1, 2, 3], acoustic wave propagation in different media [4, 5, 6], electromagnetic wave scattering as well as electrostatics [7, 8, 9, 10, 11]. An interesting aspect is that, in many cases, Robin conditions allow us to model the presence of an unresolved thin boundary layer and present themselves as a small perturbation of a Dirichlet (or Neumann) condition. This hides the asymptotic behavior of the medium being modeled within such layer, e.g., see [12] and references therein for a review. In more complex settings, e.g., as fluid-structure interaction, the Robin condition is also known to play a major role [13, 14]. An additional very appealing use of Robin conditions is their potential for managing the variability of the boundary conditions in space and time – with a single implementation – as shown in [15].

In this paper, we are interested in studying high-order unfitted methods, and in particular Robin formulations, that can provide high-order accuracy on unfitted meshes. This can be used to switch from one boundary condition to another with a single implementation, which is consistent with asymptotic approximations possibly embedded in the PDE. The focus of the paper is on the Poisson equation with a reaction term as this is the governing spatial operator in many of the physical applications mentioned above.

*High-order hp-spectral element methods.* Many problems are well known to be more efficiently approximated by high-order methods. In this paper, we target a very high-order method with convergence rates much higher than second-order. The known gains obtainable with such methods are associated with faster convergence if the solution has sufficient smoothness or appropriate *hp*-adaptive strategies are adopted [16, 17, 18].

Among high-order numerical discretizations, the spectral element method (SEM) [18], or *hp*-FEM [19], are classical choices. The SEM is a combination of two classes of methods: a high-order accurate spectral collocation method and a finite element method that is designed to handle geometry in a flexible way. It was introduced by Patera in the 1980s for solving laminar flows [20]. If the solution to the mathematical problem stated in terms of PDEs is sufficiently smooth, the errors will decay at an exponential/spectral rate when increasing the approximation degree of the polynomial basis, *p*. Also, compared to spectral collocation methods, the resulting linear system of equations has a high degree of sparsity due to the domain decomposition of the domain into elements. Ultimately, this provides efficient linear scaling opportunities when dealing with complex geometry, e.g., through geometric *p*-multigrid methods [21, 22]. The SEM also allows for the use of simplicial meshes, which have the benefit of great geometrical flexibility where necessary. In this work, SEM approximations up to an order of 11 will be considered since this is assumed to cover most practical needs.

*High order unfitted approximation and shifted boundary approach.* As already mentioned, unfitted methods avoid cumbersome mesh generation and mesh management issues in applications involving complex geometries or moving boundaries. There exist many different approaches to constructing such methods. A widely used approach is based on volume penalization, or forcing, which originates from the work of Peskin in the 1970s [23]. The method is quite robust and has been applied in many applications (see the recent review [24]) and to impose several types of boundary conditions, including Neumann and Robin [25, 26, 27]. However, its accuracy is often limited to first or second-order, which requires advanced use of adaptive techniques to offset the associated error [24, 28]. The finite element setting has shown a greater potential to allow for the construction of higher accuracy unfitted methods. Starting with the early works on GFEM, XFEM, and cut finite elements [29, 30, 31], many developments have allowed us to tackle some of the underlying initial limitations of methods based on mesh intersections. The latter may lead to geometrically ill-posed situations, which may translate into bad conditioning of the resulting problem when elliptic equations are considered. Several approaches have emerged over the years involving either the use of appropriate stabilization oper-
ators, cell aggregation, or the use of level sets to redefine the finite element spaces [32, 33, 34, 35, 36, 37, 38]. Many of these approaches allow high-degree polynomial approximations as well.

In this paper, we focus on a numerical discretization strategy for handling domain boundary conditions introduced relatively recently by G. Scovazzi and his collaborators: the shifted boundary method (SBM). The SBM was initially proposed by Main & Scovazzi in [39, 40] for solving boundary value problems. The main idea behind the SBM is first to relocate – or *shift* – the domain boundary from the *true* to an approximate – or *surrogate* – boundary. The latter is composed of appropriately defined mesh faces. The subset of mesh elements within the surrogate boundary constitutes the surrogate domain on which the discrete solution is computed. To ensure that the boundary conditions are imposed with the desired accuracy, the second idea in the SBM is to modify the problem data to account for the geometrical shift of the boundary. This step is classically used as a *shift operator* – from the surrogate to the true boundary – provided by a truncated Taylor series development. This approach provides a simple, efficient, and robust method that overcomes complex meshing tasks, naturally incorporates the ability to deal with curved geometries, and avoids potential problems with small cut cells. The method has been successfully extended to hyperbolic and parabolic equations [41, 42, 43, 44], to the management of fixed and moving interfaces [43, 45, 46, 47], to solid mechanics [48, 49], and to higher orders [42, 50, 51, 52].

*Accuracy and conditioning.* Two key elements in most of the unfitted finite element approaches mentioned are: i) the way in which the boundary conditions are enforced and ii) the polynomial approximation used to achieve such enforcement. In the context of the SBM, the data is typically extrapolated from the surrogate boundary to the true boundary. For the original formulation, this is done in terms of truncated Taylor series, which – as shown in [42] – is exactly equivalent to the evaluation of the cell polynomials onto the true boundary. This leaves several options, as the true boundary may fall within the surrogate domain or outside. The latter approach is usually the one considered in most SBM publications, with the exception of the recent work considering a mix of the two combined with adaptive mesh refinement [53]. When considering polynomials of high degree, this aspect may have a strong impact on the accuracy and stability of SBM. Extrapolation of such polynomials is known to converge less optimal and lead to less well-conditioned algebraic equations [54].

Finally, the form of the terms included in the variational statement may have an impact on the accuracy and con-
ditioning of the problem. In many of the works mentioned on high-order unfitted finite elements, including the SBM,
the boundary condition is enforced weakly. For Dirichlet conditions, this is often achieved using Nitsche's variation
formulation [55]. To impose more general conditions, generalized variants of Nitsche's method could be used, as
those proposed in [56, 57]. To our knowledge, the only example of this in the context of unfitted finite elements is
in [58]. These formulations involve several additional terms, including adjoint consistent corrections and penalization
terms. Another related and classical approach to enforce Dirichlet conditions weakly is the one proposed by Aubin,
also known as a penalty method [59]. This approach is somewhat simpler compared to Nitsche's method. However, it
is not adjoint-consistent, which may lead to sub-optimal convergence, e.g., see [60] and references therein.

1.1 Paper aim and contribution

In this paper, we propose and evaluate an SEM-SBM approach with general weak Robin boundary conditions. The first objective is – besides providing a high-order approximation of genuine Robin conditions – to obtain formulations allowing the consideration of spatially varying conditions going from Dirichlet to Neumann. Secondly, for problems involving small perturbations of one of the above limits, we would like the discrete approximation to be consistent with the asymptotics underlying the continuous PDE. These objectives will be achieved by combining the SBM with appropriate generalized weak forms [57]. Finally, we perform a thorough study of the impact of the polynomial approximation on the surrogate boundary and of the variational form used on the numerical convergence and condition number. To this end, we study three variants of the SBM involving different maps from the mesh to the surrogate boundary and generalized weak forms of both the Nitsche and Aubin types.

## 1.2 Structure of paper

The remainder of this paper is outlined as follows: In Section 2, the mathematical problem is stated for the two-
dimensional scalar Poisson problem in strong form with Dirichlet, Neumann, and Robin boundary data. Then, in
Section 3, we outline various preliminaries for the proposed model, including the computational domains, mappings,
finite element tessellation, and more. In Section 4, the numerical method – in terms of the SEM – is outlined on confor-
mal meshes. This section is followed by Section 5, where the polynomial corrected shifted boundary is introduced in
the Dirichlet and Neumann context. Hereafter, in Section 6, the case of Robin conditions are handled. In Section 7, we
present multiple numerical experiments, including convergence studies and numerical analysis of matrix conditioning.
Lastly, in Section 8, a conclusion is formed to summarize the main findings and achievements of the paper.
## 2 Mathematical problem

First, we define the mathematical problem by considering a two-dimensional ($d = 2$) planar domain in a Euclidean space, $\Omega \subset \mathbb{R}^d$, with a potentially curved Lipschitz boundary, $\Gamma = \partial\Omega$. On the boundary, we define the outward-facing normal and tangent unit vectors by $\mathbf{n} = \mathbf{n}(x, y) : \Gamma \mapsto \mathbb{R}^d$ and $\mathbf{t} = \mathbf{t}(x, y) : \Gamma \mapsto \mathbb{R}^d$, respectively. See Figure 1 for a visualization. In $\Omega$, we consider the scalar solution, $u = u(x, y) : \Omega \mapsto \mathbb{R}$, of the following elliptic problem with different non-homogeneous boundary conditions, where $u \in C^2(\Omega)$, such that

$$ 
\begin{array}{l}
-\nabla^2 u + \alpha u = f, \quad \text{in } \Omega, \\
u = u_D, \quad \text{on } \Gamma_D, \\
\nabla u \cdot \mathbf{n} = q_N, \quad \text{on } \Gamma_N, \\
u + \varepsilon(\nabla u \cdot \mathbf{n}) = u_{RD} + \varepsilon q_{RN}, \quad \text{on } \Gamma_R.
\end{array} 
\quad (1) $$

Here $\Gamma$ is partitioned into boundaries for Dirichlet ($\Gamma_D$), Neumann ($\Gamma_N$), and Robin data ($\Gamma_R$) such that $\Gamma_D \cup \Gamma_N \cup \Gamma_R = \Gamma$. Furthermore, $\alpha \ge 0$ governs the reaction term to ensure a well-posed pure Neumann problem when $\Gamma_D = \Gamma_R = \emptyset$. Also, $f = f(x, y) : \Omega \mapsto \mathbb{R}$ is a given forcing function. The Dirichlet and Neumann data are denoted by $u_D = u_D(x, y) : \Gamma_D \mapsto \mathbb{R}$ and $q_N = q_N(x, y) : \Gamma_N \mapsto \mathbb{R}$, respectively. Moreover, the Dirichlet and Neumann contributions to the Robin boundary condition are denoted $u_{RD} = u_{RD}(x, y) : \Gamma_R \mapsto \mathbb{R}$ and $q_{RN} = q_{RN}(x, y) : \Gamma_R \mapsto \mathbb{R}$, respectively. Note that, in practice, only one of the two is known. Considering both terms, however, is interesting in view of the study of the asymptotic limits.

### 2.1 Formal asymptotics of the full Robin problem

Let us focus for a moment on the Robin problem obtained when $\Gamma_R = \Gamma$ as

$$ 
\begin{array}{l}
-\nabla^2 u + \alpha u = f, \quad \text{in } \Omega, \\
u + \varepsilon(\nabla u \cdot \mathbf{n}) = u_{RD} + \varepsilon q_{RN}, \quad \text{on } \Gamma.
\end{array} 
\quad (2) $$

Interestingly, the Robin condition can be used to mimic Dirichlet and Neumann conditions as a *all-in-one* boundary formulation by choosing either small or large values of $\varepsilon$. This fact has already been exploited in [15] to impose variable conditions on conformal meshes. Here, we try to go a little further. Indeed, as recalled in the introduction, many applications involve small perturbations of a Dirichlet or Neumann condition. The cases are obtained when $\varepsilon \ll 1$ or $\varepsilon^{-1} \ll 1$ – but not zero – in (2), respectively. Using classical formal asymptotics, we can show the following.

**Proposition 2.1 (Robin-Dirichlet asymptotics).** For a smooth boundary $\Gamma$, and $u_{RD}$ and $q_{RN}$ smooth enough functions, the full Robin problem in (2) admits the following formal asymptotics for $\varepsilon \to 0$ as a series expansion as

$$ u = u_0 + \varepsilon u_1 + \sum_{m>1} \varepsilon^m u_m, \quad (3a) $$

where the leading order term, $u_0 = u_0(x, y) : \Omega \mapsto \mathbb{R}$, verifies the elliptic equation given Dirichlet boundary data as

$$ -\nabla^2 u_0 + \alpha u_0 = f, \quad \text{in } \Omega \quad \text{with} \quad u_0 = u_{D,0} := u_{RD} \quad \text{on } \Gamma, \quad (3b) $$

while the Dirichlet problems give the first and higher-order corrections $u_m = u_m(x, y) : \Omega \mapsto \mathbb{R}$ for $\forall m \ge 1$ as

$$ -\nabla^2 u_1 + \alpha u_1 = 0, \quad \text{in } \Omega \quad \text{with} \quad u_1 = u_{D,1} := q_{RN} - \partial_n u_0 \quad \text{on } \Gamma, \quad (3c) $$

and $\forall m > 1$

$$ -\nabla^2 u_m + \alpha u_m = 0, \quad \text{in } \Omega \quad \text{with} \quad u_m = u_{D,m} := -\partial_n u_{m-1} \quad \text{on } \Gamma, \quad (3d) $$

where $\partial_n$ denotes the derivative with respect to the normal direction of the boundary.

**Proposition 2.2 (Robin-Neumann asymptotic).** For a smooth boundary $\Gamma$, and $u_{RD}$ and $q_{RN}$ smooth enough functions, the full Robin problem in (2) admits the following formal asymptotics for $\varepsilon^{-1} \to 0$

$$ u = u_0 + \varepsilon^{-1} u_1 + \sum_{m>1} \varepsilon^{-m} u_m, \quad (4a) $$

where the leading order term $u_0 = u_0(x, y) : \Omega \mapsto \mathbb{R}$ verifies the elliptic equation given Neumann boundary data as

$$ -\nabla^2 u_0 + \alpha u_0 = f, \quad \text{in } \Omega \quad \text{with} \quad \partial_n u_0 = q_{N,0} := q_{RN} \quad \text{on } \Gamma, \quad (4b) $$

while the Neumann problems give the first and higher-order corrections $u_m = u_m(x, y) : \Omega \mapsto \mathbb{R}$ for $\forall m \ge 1$ as

$$ -\nabla^2 u_1 + \alpha u_1 = 0, \quad \text{in } \Omega \quad \text{with} \quad \partial_n u_1 = q_{N,1} = u_{N,1} := u_{RD} - u_0 \quad \text{on } \Gamma, \quad (4c) $$

and $\forall m > 1$

$$ -\nabla^2 u_m + \alpha u_m = 0, \quad \text{in } \Omega \quad \text{with} \quad \partial_n u_m = q_{N,m} = u_{N,m} := -u_{m-1} \quad \text{on } \Gamma. \quad (4d) $$
<figure><img src="image_3.png" /><figcaption>Figure 1: Conceptual layout of the true domain and the surrogate domain (top left). Conformal affine meshing (top right). Extrapolation mode (bottom left) and interpolation mode (bottom right). Surrogate domain boundary (solid black): $\bar{\Gamma}_h$, true domain boundary (dashed red): $\Gamma$, and triangulated mesh (dotted blue in the background): elements $\Omega^n$.</figcaption></figure>

The above propositions can be proven – for smooth enough data and domain – using classical matched asymptotic expansions in $\epsilon$ or $\epsilon^{-1}$, e.g., see [61] chapter 7 or [62] chapter 3. We omit the details for brevity. Moreover, we will say that a numerical approximation of (2) is asymptotic preserving (AP) if the discrete asymptotics of such approximation is consistent with the continuous ones at all orders in powers of $\epsilon$. Note that when the Robin condition is associated with a physical process, usually $q_{RN} = 0$, for the small $\epsilon$ limit, or $u_{RD} = 0$ in the small $\epsilon^{-1}$ limit. However, keeping these values allows us to include the discussion of when using Robin conditions to handle Dirichlet or Neumann boundaries.

## 3 Domain discretization: mesh and boundary mapping

In this section, different domain discretizations will be defined to form the geometrical understanding and notation. For this, we consider regular (quasi-uniform) tessellations using $N_{elm}$ conforming non-overlapping affine triangular elements. The computational domain, $\Omega_h$, is thus defined by

$$ \Omega_h = \bigcup_{n=1}^{N_{elm}} \Omega^n, \qquad (5) $$

with $\Omega^n$ being the $n$'th element. We denote the maximum edge length for each element as $h_{max}^n$. Similar local quantities, such as the minimum, $h_{min}^n$, and the average, $h_{avg}^n$, edge length can be computed. The associated global measure of the maximum edge length is taken as $h_{max} = \max(\{h_{max}^n\}_{n=1}^{N_{elm}})$, where $h_{min}$ and $h_{avg}$ are found similarly.

Let $\Gamma_h$ denote the boundary of the computational domain. We will consider two main cases in this paper: a conformal (fitted) and a non-conformal (unfitted) boundary approach. In the fitted case, $\Gamma_h$ is conformal with the analytical boundary $\Gamma$, and is usually defined by a classical iso-parametric high-degree finite element approximation. The analytical boundary is embedded into a non-conformal background mesh in the unfitted case. Here, we distinguish between the true domain, defined by the analytical boundary $\Gamma$, and the surrogate domain. The surrogate domain is denoted by $\bar{\Gamma}_h$. It is defined here using level set theory [63] equipped with a signed distance function – such as the Euclidean metric – to capture the geometry of the domain implicitly.
For a given element, $\Omega^n$, of arbitrary type, three situations can occur depending on the level set sign at its spanning vertices: inside $\Omega$ (all positive signs), outside $\Omega$ (all negative signs), crossed by the true boundary $\Gamma$ (mix of signs). Once the elements crossed by the true boundary are known, the set of element edges closest to the true boundary can be located using a closest-point map (or some other projection strategy). Ultimately, we obtain a set of edges defining the surrogate boundary, $\bar{\Gamma}_h$. This implicitly defines a surrogate domain $\bar{\Omega}_h$ composed of the elements contained within the surrogate boundary. As before, we denote these elements by $\{\Omega^n\}_{n=1}^{N_{\text{elm}}}$. On $\bar{\Gamma}_h$ we define the outward-facing normal unit vector by $\bar{n} = \bar{n}(x, y) : \bar{\Gamma}_h \to \mathbb{R}^d$.

## 3.1 Boundary mapping

For the purposes of this paper, we consider two different ways of defining the surrogate boundary $\bar{\Gamma}_h$. These are referred to as extrapolation and interpolation mode. These are defined as follows

* Extrapolation: $\bar{\Omega}_h$ is constructed using only whole elements inside the true domain $\Omega$.

* Interpolation: $\bar{\Omega}_h$ is constructed by including all whole elements inside $\Omega$ and the ones intersecting $\Gamma$.

To support the definitions, the true and different surrogate domains are illustrated in a Cartesian coordinate system spanned by the horizontal $x$-axis and the vertical $y$-axis in Figure 1. Also, an example of a conformal affine mesh is visualized in the figure. Interpolation-based embedded methods have been studied in recent works, see, e.g., [53, 54], as well as [53] in the context of the second-order SBM approximations on adaptive quadrilateral meshes. In this work, we consider the latter approach, combined with very high-order spectral element approximations on simplicial meshes, although the exact same formulations can be used on tensor product Cartesian meshes.

**Remark 3.1.** In the interpolation mode, the discrete domain includes a small region outside the analytical one. This requires the definition of an extension of the problem data and possibly of the exact solution in a small region outside the physical domain. This is a well-known issue related to the construction of embedded methods (see, e.g., [37, 64, 65, 66, 67]). Some practical solutions have been proposed in the past, e.g., in [64, 68]. We do not dwell on this issue in this work. All the applications considered involve manufactured solutions, for which this extension is naturally obtained from the (known) analytical solution.

Once the surrogate boundary is defined, a key ingredient is the availability of a mapping from the surrogate to the true boundary, where the boundary conditions are defined. This mapping is denoted by $\mathcal{M}$ as

$$ \mathcal{M}(\bar{x}) : \bar{\Gamma}_h \mapsto \Gamma, \quad \bar{x} \mapsto x, \qquad (6) $$

where a point on $\Gamma$ is denoted by $x = (x, y)$ and similarly on $\bar{\Gamma}_h$ by $\bar{x} = (x, y)$. Using the above, we have that $x = \mathcal{M}(\bar{x})$, i.e., the set of points on the true boundary spanned by the points of $\bar{\Gamma}_h$. This allows us to define a distance vector as $d = d(\bar{x}) = x - \bar{x}$. Note that there are several possible ways to define the mapping. For example, with the closest-point projection, any point on $\bar{\Gamma}$ is mapped to $\Gamma$ using the true normal, $n$, such that $d \cdot n \equiv 1$ and hereby $d = \|d\|n$. Other definitions are explored in this paper. The mapping, $\mathcal{M}$, allows us to define extensions of any particular function from the surrogate to the true boundary (and vice-versa).

# 4 Spectral element discretization on conformal meshes

In the following, we consider and define the spectral element method on conformal meshes for generalized Aubin's and Nitsche's methods. This is then related to the asymptotic preserving properties of the two variational forms.

## 4.1 Spectral element expansion

We consider approximate solutions in the finite element space of continuous and piece-wise polynomials of degree $P$ as

$$ \mathcal{V}_h^P = \{v \in C^0(\Omega); \forall n \in \{1, \dots, N_{\text{elm}}\}, v|_{\Omega^n} \in \mathbb{P}^P\}. \qquad (7) $$

Thus, the discrete unknown, $u_h$, can be expressed as $N_{\text{elm}}$ local element-based solutions, $u_h^n = u_h^n(x, y) : \Omega^n \mapsto \mathbb{R}$, as

$$ u_h = u_h(x, y) = \bigoplus_{n=1}^{N_{\text{elm}}} u_h^n(x, y), \quad \text{for } (x, y) \in \Omega^n, \qquad (8) $$
<figure><img src="image_4.png" /><figcaption>Figure 2: Concept of the local element domain, $\Omega^n$, the reference domain, $\mathcal{R}$, and the affine mapping, $\Psi: \Omega^n \to \mathcal{R}$.</figcaption></figure>

In this work, $u_h$ is continuous and expanded on local bases whose restriction to an element $\Omega^n$ may be equivalently defined in terms of nodal unknowns or modal expansions as in [69]

$$u_h^n(x, y) \approx \sum_{m=1}^{N_{ep}} \hat{u}_m^n \psi_m(x, y) = \sum_{i=1}^{N_{ep}} u_i^n h_i(x, y), \quad \text{for } (x, y) \in \Omega^n, \qquad (9)$$

where $N_{ep}$ is the number of quadrature element points in an element. The local modal form is expressed as a sum of $N_{ep}$ local expansion coefficients, $\hat{u}_m^n$, and associated modal two-dimensional basis functions, $\psi_m$, from the family of Jacobi polynomials up to order $P$ defined on the $N_{ep}$ quadrature points. Moreover, the local nodal counterpart is expressed as a sum of $N_{ep}$ local solution values, $u_i^n$, and nodal two-dimensional basis functions, $h_i$, as Lagrange polynomials up to order $P$ defined on the $N_{ep}$ quadrature points. For completeness, we denote the $N_{ep}$ quadrature points in $\Omega^n$ as $\{x_i, y_i\}_{i=1}^{N_{ep}}$, such that $u_i^n = u^n(x_i, y_i)$ and the Lagrange polynomials have the cardinal property $h_j(x_i, y_i) = \delta_{ij}$. In the case of finite triangular elements, the number of element points, $N_{ep}$, is linked to the polynomial order of the basis functions, $P$, via $N_{ep} = \frac{P+1}{2}$.

In practice, we adopt an iso-parametric approach, where each physical element, $\Omega^n$, is mapped onto a triangular reference domain, $\mathcal{R} = \{\mathbf{r} = (r, s) | (r, s) \ge -1; r + s \le 0\}$, by introducing the affine mapping, $\Psi = \Psi(r, s): \Omega^n \to \mathcal{R}$, such that $(x, y) = \Psi(r, s)$ and $(r, s) = \Psi^{-1}(x, y)$. See Figure 2 for a visualization. In accordance with [70], we construct an orthonormal basis, $\psi_m(r, s)$, on a set of quadrature points, $(r, s)$, in $\mathcal{R}$. The basis functions are chosen as the family of Jacobi polynomials. Moreover, the quadrature points are chosen to be Gauss-Lobatto distributed by applying an optimized blending and wrapping procedure of the nodes by [71] such that the Lebesque constant is minimized. This delicate arrangement ultimately leads to a well-behaved generalized Vandermonde matrix for interpolation as

$$\mathcal{V} = \mathcal{V}_{ij} = \psi_j(r_i, s_i), \quad \text{for } (i, j) = 1, \dots, N_{ep}, \quad \mathcal{V} \in \mathbb{R}^{N_{ep} \times N_{ep}}. \qquad (10)$$

## 4.2 Variational formulation on conformal meshes: a generalized Nitsche's method

We consider a weak form of (1) based on a generalized Nitsche's method, hereby embedding all types of boundary conditions. On conforming meshes, the method's consistency, coercivity, and error estimates have been widely studied in the past, e.g., see [55, 56, 57]. The classical $L^2$ scalar product over the set $A$ is denoted by

$$(f, g)_A = \int_A f g \, dA. \qquad (11)$$

The weak formulation is now to find $u_h \in \mathcal{V}_h^P$ such that $\forall v_h \in \mathcal{V}_h^P$, where the following holds

$$\begin{align} (\nabla u_h, \nabla v_h)_{\Omega_h} &+ \alpha(u_h, v_h)_{\Omega_h} - (f, v_h)_{\Omega_h} - (q_N, v_h)_{\Gamma_{h,N}} \nonumber \\ & - (\nabla u_h \cdot \mathbf{n}, v_h)_{\Gamma_{h,D}} - (u_h - u_D, \nabla v_h \cdot \mathbf{n})_{\Gamma_{h,D}} + \gamma^{-1}(u_h - u_D, v_h)_{\Gamma_{h,D}} \nonumber \\ & - \frac{\gamma}{\gamma + \varepsilon} (\nabla u_h \cdot \mathbf{n}, v_h)_{\Gamma_{h,R}} - \frac{\gamma}{\gamma + \varepsilon} (u_h - u_{RD}, \nabla v_h \cdot \mathbf{n})_{\Gamma_{h,R}} + \frac{1}{\gamma + \varepsilon} (u_h - u_{RD}, v_h)_{\Gamma_{h,R}} \tag{12} \\ & - \frac{\varepsilon}{\gamma + \varepsilon} (q_{RN}, v_h)_{\Gamma_{h,R}} - \frac{\varepsilon\gamma}{\gamma + \varepsilon} (\nabla u_h \cdot \mathbf{n} - q_{RN}, \nabla v_h \cdot \mathbf{n})_{\Gamma_{h,R}} = 0. \nonumber \end{align}$$

The variational form above, which includes several stabilizing and symmetrizing terms, has been analyzed in depth in the past and shown to be continuous and coercive [55, 56, 57]. In particular, the parameter $\gamma^{-1}$ is the penalty
coefficient of the classical Dirichlet Nitsche method and of the generalized Nitsche method of [57]. These have been studied in several past works, e.g., see [56, 57, 72] and references therein.

With some algebra, the Euler-Lagrange equations – and its strong consistency – can be shown to be

$$(-\nabla^2 u_h + \alpha u_h - f, v_h)_{\Omega_h} + (u_h - u_D, \gamma^{-1}v_h - \nabla v_h \cdot \mathbf{n})_{\Gamma_{h,D}} + (\nabla u_h \cdot \mathbf{n} - q_N, v_h)_{\Gamma_{h,N}} + \frac{1}{\gamma + \varepsilon}(u_h - u_{RD} + \varepsilon\nabla u_h \cdot \mathbf{n} - \varepsilon q_{RN}, v_h - \gamma\nabla v_h \cdot \mathbf{n})_{\Gamma_{h,R}} = 0. \quad (13)$$

To simplify the notation, from now on, we will omit the subscript "h", and assume implicitly that all quantities involving discrete polynomial functions are expressed and integrated on mesh entities. Also, we recast the variational form (12) in the short notation as

$$a(u, v) = b(v), \quad (14a)$$

where

$$a(u, v) = a(u, v)_{\Omega} + a(u, v)_{\Gamma_D} + a(u, v)_{\Gamma_N} + a(u, v)_{\Gamma_R}, \quad (14b)$$

$$b(v) = b(v)_{\Omega} + b(v)_{\Gamma_D} + b(v)_{\Gamma_N} + b(v)_{\Gamma_R}, \quad (14c)$$

with

$$a(u, v)_{\Omega} = (\nabla u, \nabla v)_{\Omega} + \alpha(u, v)_{\Omega}, \quad (14d)$$

$$a(u, v)_{\Gamma_D} = \gamma^{-1}(u, v)_{\Gamma_D} - (\nabla u \cdot \mathbf{n}, v)_{\Gamma_D} - (u, \nabla v \cdot \mathbf{n})_{\Gamma_D}, \quad (14e)$$

$$a(u, v)_{\Gamma_N} = 0, \quad (14f)$$

$$a(u, v)_{\Gamma_R} = \frac{1}{\varepsilon + \gamma}(u, v)_{\Gamma_R} - \frac{\gamma}{\varepsilon + \gamma}(\nabla u \cdot \mathbf{n}, v)_{\Gamma_R} - \frac{\gamma}{\varepsilon + \gamma}(u, \nabla v \cdot \mathbf{n})_{\Gamma_R} - \frac{\varepsilon\gamma}{\varepsilon + \gamma}(\nabla u \cdot \mathbf{n}, \nabla v \cdot \mathbf{n})_{\Gamma_R}, \quad (14g)$$

$$b(v)_{\Omega} = (f, v)_{\Omega}, \quad (14h)$$

$$b(v)_{\Gamma_D} = \gamma^{-1}(u_D, v)_{\Gamma_D} - (u_D, \nabla v \cdot \mathbf{n})_{\Gamma_D}, \quad (14i)$$

$$b(v)_{\Gamma_N} = (q_N, v)_{\Gamma_N}, \quad (14j)$$

$$b(v)_{\Gamma_R} = \frac{1}{\varepsilon + \gamma}(u_{RD}, v)_{\Gamma_R} + \frac{\varepsilon}{\varepsilon + \gamma}(q_{RN}, v)_{\Gamma_R} - \frac{\gamma}{\varepsilon + \gamma}(u_{RD}, \nabla v \cdot \mathbf{n})_{\Gamma_R} - \frac{\varepsilon\gamma}{\varepsilon + \gamma}(q_{RN}, \nabla v \cdot \mathbf{n})_{\Gamma_R}. \quad (14k)$$

## 4.3 Variational formulation on conformal meshes: a generalized Aubin's method

We now consider a reduced form of Nitche's formulation essentially generalizing the method by Aubin [59]. By inspecting (13), we can see that the Robin contributions can be written as a linear combination of a Neumann penalty with coefficient $\varepsilon/(\varepsilon+\gamma)$, a Nitche penalty with penalty parameter $\gamma^{-1}$, and coefficient $\gamma/(\varepsilon+\gamma)$. To obtain a method similar to Aubin's penalty formulation, we consider the same above combination without the symmetrizing terms in the penalties. In the short-form notation, this gives the generalized Aubin variational form as

$$a(u, v)_{\Omega} = (\nabla u, \nabla v)_{\Omega} + \alpha(u, v)_{\Omega}, \quad (15a)$$

$$a(u, v)_{\Gamma_D} = \gamma^{-1}(u, v)_{\Gamma_D} - (\nabla u \cdot \mathbf{n}, v)_{\Gamma_D}, \quad (15b)$$

$$a(u, v)_{\Gamma_N} = 0, \quad (15c)$$

$$a(u, v)_{\Gamma_R} = \frac{1}{\varepsilon + \gamma}(u, v)_{\Gamma_R} - \frac{\gamma}{\varepsilon + \gamma}(\nabla u \cdot \mathbf{n}, v)_{\Gamma_R}, \quad (15d)$$

and

$$b(v)_{\Omega} = (f, v)_{\Omega}, \quad (15e)$$

$$b(v)_{\Gamma_D} = \gamma^{-1}(u_D, v)_{\Gamma_D}, \quad (15f)$$

$$b(v)_{\Gamma_N} = (q_N, v)_{\Gamma_N}, \quad (15g)$$

$$b(v)_{\Gamma_R} = \frac{1}{\varepsilon + \gamma}(u_{RD}, v)_{\Gamma_R} + \frac{\varepsilon}{\varepsilon + \gamma}(q_{RN}, v)_{\Gamma_R}. \quad (15h)$$

Integration by part yields the Euler-Lagrange equations – that emphasizes the consistency of the formulation – as

$$(-\nabla^2 u + \alpha u - f, v)_{\Omega} + \gamma^{-1}(u - u_D, v)_{\Gamma_D} + (\nabla u \cdot \mathbf{n} - q_N, v)_{\Gamma_N} + \frac{1}{\varepsilon + \gamma}(u - u_{RD} + \varepsilon \nabla u \cdot \mathbf{n} - \varepsilon q_{RN}, v)_{\Gamma_R} = 0. \quad (16)$$
Note that the coercivity of (15a)-(15d) can be shown by classical techniques based on trace inequalities to bound the gradient term on the boundary, e.g., see [59, 73] and [74, Chapter 37]. As mentioned in other works, in the Dirichlet case, e.g., see [60], the penalty method requires appropriate choices of the penalty constant, and is in general less consistent with respect to the symmetric and adjoint consistent Nitsche method. However, as we will see in Section 7, very high orders of approximation on embedded boundaries and the presence of additional derivative terms have a strong impact on the conditioning, which makes this weak form of interest to consider in detail.

## 4.4 Algebraic system of equations

Finally, following the standard Galerkin spectral element approach, a system of equations is obtained by choosing the test function, $v$, to equal the nodal basis functions, $h$, outlined in Section 4.1. Hereafter, we evaluate all the integrals using accurate numerical quadrature [18, 20, 69]. This ultimately leads to a large algebraic system of equations for the nodal unknowns as

$$Au = b, \quad \text{where} \quad A \in \mathbb{R}^{N_{\text{DOF}} \times N_{\text{DOF}}}, \quad \text{and} \quad (u, b) \in \mathbb{R}^{N_{\text{DOF}}}. \qquad (17)$$

with $N_{\text{DOF}}$ being the number of degrees of freedom (DOF) – corresponding to the number of nodal points, expansion coefficients, and/or basis functions – in the discrete scheme. The system matrix and vector are denoted by $A$ and $b$, respectively, corresponding to the weak statements’ left- and right-hand sides. The discrete solution vector $u$ has been obtained using a direct linear solver to invert $A$. Contributions to $A$ and $b$ are obtained on a local element level with discrete operators for integration and differentiation as defined in [70].

## 4.5 Asymptotic preserving properties of the Robin variational forms

The boundary entry of most interest in (12) is the one related to the Robin terms. Because of their asymptotic properties, described by Proposition 2.1 and 2.2, these terms allow setting – with a single implementation – different types of boundary conditions by just choosing a different value of the $\varepsilon$ parameter. To the authors’ knowledge, this approach has never been tested in the context of embedded methods and – in any case – not in the framework of the SBM.

For this to be possible, it is important that the discrete equations correctly mimic the asymptotic behavior of the continuous problem. In this case, we speak of asymptotic preserving approximations. It has already been remarked in [57] that the classical implementation of the Robin conditions (based on a modified Neumann flux) does not lend itself to this purpose. They show numerically that their proposed method is well posed in both the zeroth of $\varepsilon$ and $\varepsilon^{-1}$. Here, we go further and prove the asymptotic consistency of the method with the problems defined by Proposition 2.1 and 2.2.

### 4.5.1 The Nitsche-Robin method

First, recall that we solely consider the Robin boundary, such that $\Gamma = \Gamma_R$. Then, we introduce the following notation

$$\mathcal{I}_{\Omega}(u) = (\nabla u, \nabla v)_{\Omega} + \alpha(u, v)_{\Omega}, \quad \mathcal{N}_{\Gamma}(u) = (u, v/\gamma - \nabla v \cdot \mathbf{n})_{\Gamma}. \qquad (18)$$

Now, we consider the generalized Nitsche’s method, which – after some manipulations – can be written as

$$\mathcal{I}_{\Omega}(u) - (f, v)_{\Omega} - (\nabla u \cdot \mathbf{n}, v)_{\Gamma} + \frac{\gamma}{\gamma + \varepsilon} \mathcal{N}_{\Gamma}(u - u_{RD}) + \frac{\varepsilon\gamma}{\gamma + \varepsilon} \mathcal{N}_{\Gamma}(\nabla u \cdot \mathbf{n} - q_{RN}) = 0. \qquad (19)$$

We will now consider the asymptotic limits, $\varepsilon \to 0$ and $\varepsilon^{-1} \to 0$, and the corresponding asymptotic developments in (3a) and (4a). We can prove the following two results in the upcoming propositions.

**Proposition 4.1** (AP property of the Nitsche-Robin method in the Dirichlet limit $\varepsilon \to 0$). The asymptotic approximations obtained when $\varepsilon \to 0$ from the full Robin problem in (19) in the finite element space $V_h^P$ for the asymptotic expansion modes $u_m$ in (3a) define a cascade of variational problems which are all coercive and consistent with (3b) through (3d). In particular, the $m$'th order can be obtained from the previous ones as the solution of

$$\mathcal{I}_{\Omega}(u_m) - (f_m, v)_{\Omega} - (\nabla u_m \cdot \mathbf{n}, v)_{\Gamma} + \mathcal{N}_{\Gamma}(u_m - u_{D,m}) = \mathcal{E}_m, \qquad (20)$$

where the left-hand side defines the classical Nitsche-Dirichlet variational form with $f_0 = f$ and $f_m = 0$ for $m \ge 1$, while the right-hand side introduces asymptotic error terms with $\mathcal{E}_0 = 0$ and

$$\mathcal{E}_m = - \sum_{\ell=0}^{m-1} (-\gamma)^{-(m-\ell)} \mathcal{N}_{\Gamma}(u_{\ell} - u_{D,\ell}), \quad \text{for } m \ge 1. \qquad (21)$$
**Proof.** Equation (20) and (21) are readily obtained by introducing the ansatz from (3a) in (19), multiplying by $\gamma + \varepsilon$, and equating – term by term – all the coefficients of $\varepsilon^m$ powers. All problems are verified by the exact modes of (3b)-(3d), which can be shown classically by integrating the domain equations by parts and using the boundary conditions. The coercivity of the classical Nitsche’s variational form has been studied widely, e.g., see [56] and references therein. □

Note that in the above development, the second Nitsche term in (19) involving the normal derivatives plays an important role when matching the correct Dirichlet values at all orders. This term affects the form of the Neumann asymptotic limit. In particular, we have the following equivalent result in the latter case.

**Proposition 4.2** (AP property of the Nitsche-Robin method in the Neumann limit $\varepsilon^{-1} \to 0$). The asymptotic approximations obtained when $\varepsilon^{-1} \to 0$ from the full Robin problem in (19) in the finite element space $V_h^P$ for the asymptotic expansion modes $u_m$ in (4a) define a cascade of variational problems which are all coercive and consistent with (4b) through (4d). In particular, the m'th order can be obtained from the previous ones as the solution of

$$ \mathcal{I}_{\Omega}(u_m) - (f_m, v)_{\Omega} - (q_{N,m}, v)_{\Gamma} - \gamma(\nabla u_m \cdot \mathbf{n} - q_{N,m}, \nabla v \cdot \mathbf{n}) = \mathcal{E}_m, \quad (22) $$

where the left-hand side defines a Neumann variational form with an additional penalty with $f_0 = f$ and $f_m = 0$ for $m \ge 1$, while the right-hand side introduces asymptotic error terms as $\mathcal{E}_0 = 0$ and

$$ \mathcal{E}_m = \sum_{\ell=0}^{m-1} (-\gamma)^{(m+1-\ell)} \mathcal{N}_{\Gamma}(\nabla u_{\ell} \cdot \mathbf{n} - q_{N,\ell}), \quad \text{for } m \ge 1. \quad (23) $$

**Proof.** Equation (22) and (23) are readily obtained by introducing ansatz from (4a) in (19), multiplying by $\varepsilon^{-1}\gamma + 1$, and equating – term by term – all the coefficients of $\varepsilon^{-m}$ powers. Integrating by parts of (22) shows immediately the exact consistency with (4b) to (4d), for all mode equations. Concerning coercivity, we make use of the trace inequality

$$ \|h^{1/2} \nabla u\|_{0,\Gamma} \le C_{\Gamma} \|\nabla u\|_{0,\Omega}. \quad (24) $$

This allows us to argue that

$$ (\nabla u_m, \nabla u_m)_{\Omega} + \alpha(u_m, u_m)_{\Omega} - \gamma(\nabla u_m \cdot \mathbf{n}, \nabla u_m \cdot \mathbf{n})_{\Gamma} \ge (1 - h^{-1}C_{\Gamma}\gamma)(\nabla u_m, \nabla u_m)_{\Omega} + \alpha(u_m, u_m)_{\Omega}, \quad (25) $$

which shows the coercivity of the variational form as soon as

$$ 1 - h^{-1}C_{\Gamma}\gamma \ge C > 0. \quad (26) $$

□

### 4.5.2 The Aubin-Robin method

We proceed as in the previous section. With a similar notation, we recast the Aubin-Robin method as

$$ \mathcal{I}_{\Omega}(u) - (f, v)_{\Omega} - (\nabla u \cdot \mathbf{n}, v)_{\Gamma} + \frac{\gamma}{\gamma + \varepsilon} \mathcal{A}_{\Gamma}(u - u_{RD}) + \frac{\varepsilon\gamma}{\gamma + \varepsilon} \mathcal{A}_{\Gamma}(\nabla u \cdot \mathbf{n} - q_{RN}) = 0, \quad (27) $$

where $\mathcal{I}_{\Omega}$ is defined as in (18), while

$$ \mathcal{A}_{\Gamma}(u) = (u, v/\gamma)_{\Gamma}. \quad (28) $$

The scheme is formally identical to the Nitsche one, except for the definition of the penalty which does not include derivatives of the test function. Proceeding as before, we can prove the following two results.

**Proposition 4.3** (AP property of the Aubin-Robin method in the Dirichlet limit $\varepsilon \to 0$). The asymptotic approximations obtained when $\varepsilon \to 0$ from the full Robin problem in (27) in the finite element space $V_h^P$ for the asymptotic expansion modes $u_m$ in (3a) define a cascade of variational problems which are all coercive and consistent with (3b) through (3d). In particular, the m'th order can be obtained from the previous ones as the solution of

$$ \mathcal{I}_{\Omega}(u_m) - (f_m, v)_{\Omega} - (\nabla u_m \cdot \mathbf{n}, v)_{\Gamma} + \mathcal{A}_{\Gamma}(u_m - u_{D,m}) = \mathcal{E}_m, \quad (29) $$

where the left-hand side defines the classical Aubin-Dirichlet variational form with $f_0 = f$ and $f_m = 0$ for $m \ge 1$, while the right-hand side introduces asymptotic error terms with $\mathcal{E}_0 = 0$ and

$$ \mathcal{E}_m = - \sum_{\ell=0}^{m-1} (-\gamma)^{-(m-\ell)} \mathcal{A}_{\Gamma}(u_{\ell} - u_{D,\ell}), \quad \text{for } m \ge 1. \quad (30) $$
*Proof.* The proof of the consistency is identical to the one of Proposition 4.1 (not reported for brevity). The coercivity part of the proof is identical to the classical coercivity proof for the Aubin-Dirichlet method. □

Similarly, in the Neumann asymptotic limit, we have the following equivalent result.

**Proposition 4.4** (AP property of the Aubin-Robin method in the Neumann limit $\varepsilon^{-1} \to 0$). The asymptotic approximations obtained when $\varepsilon^{-1} \to 0$ from the full Robin problem in (27) in the finite element space $V_h^P$ for the asymptotic expansion modes $u_m$ in (4a) define a cascade of variational problems which are all coercive and consistent with (4b) through (4d). In particular, the m'th order can be obtained from the previous ones as the solution of

$$ \mathcal{I}_{\Omega}(u_m) - (f_m, v)_{\Omega} - (q_{N,m}, v)_{\Gamma} = \mathcal{E}_m, \quad (31) $$

where the left-hand side defines a Neumann variational form with an additional penalty with $f_0 = f$ and $f_m = 0$ for $m \ge 1$, while the right-hand side introduces asymptotic error terms as $\mathcal{E}_0 = 0$ and

$$ \mathcal{E}_m = \sum_{\ell=0}^{m-1} (-\gamma)^{(m+1-\ell)} A_{\Gamma}(\nabla u_{\ell} \cdot \mathbf{n} - q_{N,\ell}), \quad \text{for } m \ge 1. \quad (32) $$

*Proof.* The proof is identical to the one of Proposition 4.2 (not reported for brevity). The coercivity part of the proof is identical to the standard Neumann method. □

### 4.5.3 On the stabilization parameter

To conclude this section, we stress the impact of the stabilization parameter $\gamma$ on the stability and consistency order. In particular, the last proposition shows that for the Neumann limit to be coercive, the natural choice would be to set

$$ \gamma = C_{\gamma}h, \quad C_{\gamma}C_{I} \le 1 - C, \quad (33) $$

for some strictly positive $C > 0$. This is somewhat in line with the fact that $\gamma^{-1}$ plays the role of the usual stabilization parameter in the Nitsche terms, $N_{\Gamma}$, and (33) corresponds to the classical $h^{-1}$ scaling for this parameter, see [56] and references therein.

As a side remark, the above choice – which is standard for Nitsche’s method – affects the magnitude of the error terms in Proposition 4.1. In particular, note that (21) involves scaling constants of the order $\gamma^{-m}$ which reduce the formal consistency with the higher order asymptotics by a factor of $h^{-m}$. This is of course only relevant when one is interested in reproducing higher ($\ge 2$) asymptotic modes.

## 5 Embedded discrete equations via shifted boundary polynomial corrections

In the following, we switch to a non-conformal mesh, whereupon we will solve (1) numerically. As discussed in Section 3.1, we will thus work on a surrogate domain, $\bar{\Omega}$, with a surrogate boundary $\bar{\Gamma}$, split into Dirichlet, Neumann, and Robin boundaries, $\bar{\Gamma}_D$, $\bar{\Gamma}_N$, and $\bar{\Gamma}_R$, respectively (recall that the subscript $h$ is omitted for simplicity). The definition of the appropriate boundary conditions is obtained via the mapping in (6). We thus define a weak form of the problem on the surrogate domain which can be abstractly written as

$$ a\langle u, v \rangle = b\langle v \rangle, \quad (34a) $$

where the bi-linear and linear operators can be expanded, respectively, as

$$ a\langle u, v \rangle = a\langle u, v \rangle_{\bar{\Omega}} + a\langle u, v \rangle_{\bar{\Gamma}_D} + a\langle u, v \rangle_{\bar{\Gamma}_N} + a\langle u, v \rangle_{\bar{\Gamma}_R}, \quad (34b) $$

$$ b\langle v \rangle = b\langle v \rangle_{\bar{\Omega}} + b\langle v \rangle_{\bar{\Gamma}_D} + b\langle v \rangle_{\bar{\Gamma}_N} + b\langle v \rangle_{\bar{\Gamma}_R}, \quad (34c) $$

where

$$ a\langle u, v \rangle_{\bar{\Omega}} = (\nabla u, \nabla v)_{\bar{\Omega}} + \alpha(u, v)_{\bar{\Omega}}, \quad (34d) $$

$$ b\langle v \rangle_{\bar{\Omega}} = (f, v)_{\bar{\Omega}}. \quad (34e) $$

To close the problem above, we need to specify boundary data on $\bar{\Gamma}$ as

$$ \bar{u}_D = u(\bar{x}), \quad \text{for } \bar{x} \in \bar{\Gamma}_D, \quad (34f) $$

$$ \bar{q}_N = \nabla u(\bar{x}) \cdot \bar{n}, \quad \text{for } \bar{x} \in \bar{\Gamma}_N, \quad (34g) $$

$$ \bar{u}_{RD} = u(\bar{x}), \quad \text{for } \bar{x} \in \bar{\Gamma}_R, \quad (34h) $$

$$ \bar{q}_{RN} = \nabla u(\bar{x}) \cdot \bar{n}, \quad \text{for } \bar{x} \in \bar{\Gamma}_R. \quad (34i) $$

We are now left with the challenge of how to define the surrogate boundary data above, given the data assigned on the true boundary $\Gamma$: $(u_D, q_N, u_{RD}, q_{RN})$ in (1). Using these values with no modification would result in a $O(\|d\|)$ error, thus limiting the method to first order.
## 5.1 Shifted boundary polynomial corrections

As discussed in the introduction, several approaches exist to correct this error between the two boundaries. In this work, we follow the shifted boundary approach introduced by Main & Scovazzi for elliptic and parabolic problems [39, 40], and later extended to various other applications [41, 45, 46, 47, 48, 49, 75]. In the original SBM formulation, a second-order accurate truncated Taylor expansion was introduced to match the data on boundaries. This approach has been shown to have a lot of potential in providing high-order approximations in embedded domains, avoiding classical drawbacks such as small cut-cell problems encountered in some cut-FEM methods or consistency limitations in immersed boundaries. The SBM only requires the knowledge of the mapping in (6) to enforce the given conditions while preserving the accuracy of the underlying finite element discretization. Some extensions beyond the second order have been proposed, initially in [50], and more recently in [42, 51, 52]. The strategy to correct the boundary data is the same in the high-order case: use a Taylor series development truncated at the proper order. As the accuracy increases, the number of derivative terms to be included in the series becomes quite large and cumbersome, especially in three spatial dimensions. A simplified formulation of the same method, avoiding the evaluation of all the Taylor basis terms at quadrature points, has been introduced in [42]. We use the same one here and extend it to the case of Neumann and Robin boundary conditions. Despite being equivalent to the Taylor series-based SBM, the formulation proposed in [42] involves a simple way to perform extrapolation using the underlying finite element approximation. In the following sections, we elaborate on the details of this technique.

### 5.1.1 Polynomial correction on the Dirichlet boundary

Following the SBM ideas, we start by matching the Dirichlet data, $u_D(x)$, with a k-th order multivariate Taylor series expansion of $u(\bar{x})$ in a neighborhood of $\bar{x} \in \bar{\Gamma}$. We consider at first the Nitsche's formulation. The adaptation of the same ideas to the penalty method is recalled later. In multinomial notation, this leads to

$$u_D(x) = u(\bar{x}) + D_D^k(d) + \mathcal{O}(\|d(\bar{x})\|^{k+1}), \quad \text{where} \quad D_D^k = \sum_{\alpha=1}^k \frac{1}{\alpha!} (\nabla^\alpha) (u) (\bar{x}) d(\bar{x})^\alpha, \quad (35)$$

where we recall that $d = x - \bar{x}$ and $\mathcal{O}(\|d(\bar{x})\|^{k+1})$ is the error associated with the truncated series. Neglecting this error, the above relation provides a high-order definition of the Dirichlet data on the surrogate boundary as

$$\bar{u}_D = u(\bar{x}) = u_D(x) - D_D^k. \quad (36)$$

This is the original SBM idea studied in, e.g., [51, 52] in the high-order case. In particular, as we will shortly see that (36) can be inserted directly in the standard Nitsche formulation to obtain a modified variational form, which, on unfitted meshes, is $(k+1)$'th order accurate. That is, without the need to compute any mesh intersections or define any cut-cells. Thus avoiding many of the issues associated with these cells. Although not strictly necessary, symmetrizing terms can also be added to improve the formal stability of the method [39, 50, 52]. In this work, we use a different view of the same method, which follows [42]. The key idea introduced in the reference is to recognize that the expression

$$p(x) = u(\bar{x}) + D_D^k(x - \bar{x}), \quad (37)$$

is simply the expansion of the polynomial $p$ onto the local Taylor basis $(1, x - \bar{x}, y - \bar{y}, (x - \bar{x})^2/2, (y - \bar{y})^2/2, (x - \bar{x})(y - \bar{y})$, etc.). The coefficients in this expansion are the derivatives of $u$, which are recovered from the underlying finite element polynomial in the standard SBM implementation, evaluated at $\bar{x}$. The point-wise value obtained with this expansion is clearly the same as the one provided by any other basis. With this, there is no reason to change the basis at each quadrature point, and one can use all along the underlying finite element polynomial expansion used in the volume terms. In particular, for a polynomial of degree $k$, the Taylor term can be simply expressed as

$$D_D^k = u(x) - u(\bar{x}). \quad (38)$$

This provides the following expression for the corrected value for the Dirichlet data on the surrogate boundary

$$\bar{u}_D = u(\bar{x}) + u_D(x) - u(x). \quad (39)$$

So, the SBM correction to match the boundary data can be written as the difference between the finite element cell polynomial evaluated on the true boundary and the data itself. The use of the available polynomial expansions avoids the cumbersome re-evaluations of all high-order derivatives. Following the original SBM idea, we now consider the standard Nitsche Dirichlet variational form with the modified data (39), which – using the notation $x = \mathcal{M}(\bar{x})$ – gives

$$\begin{align}
a\langle u, v \rangle_{\bar{\Gamma}_D} &= \gamma^{-1}(u(x), v)_{\bar{\Gamma}_D} - (\nabla u \cdot \bar{n}, v)_{\bar{\Gamma}_D} - (u(x), \nabla v \cdot \bar{n})_{\bar{\Gamma}_D}, \tag{40} \\
b\langle v \rangle_{\bar{\Gamma}_D} &= \gamma^{-1}(u_D, v)_{\bar{\Gamma}_D} - (u_D, \nabla v \cdot \bar{n})_{\bar{\Gamma}_D}. \nonumber
\end{align}$$
Note that if $u(x)$ is replaced by $u(\bar{x}) + D_D^k$ with $D_D^k$ given by (35), we obtain the usual SBM formulation, e.g., see [51]. However, additional terms are added to enhance the symmetry of $a\langle u, v\rangle_{\bar{\Gamma}_D}$. In the classical SBM formulation, these terms involve the Taylor series terms applied to the test function. The symmetrized form can be written as

$$ 
\begin{aligned}
a\langle u, v\rangle_{\bar{\Gamma}_D} &= \gamma^{-1}(u(x), v(x))_{\bar{\Gamma}_D} - (\nabla u \cdot \bar{n}, v)_{\bar{\Gamma}_D} - (u(x), \nabla v \cdot \bar{n})_{\bar{\Gamma}_D}, \\
b\langle v\rangle_{\bar{\Gamma}_D} &= \gamma^{-1}(u_D, v(x))_{\bar{\Gamma}_D} - (u_D, \nabla v \cdot \bar{n})_{\bar{\Gamma}_D}.
\end{aligned} 
\quad (41) $$

Note that, for a full Dirichlet problem, the resulting method is exactly the same as the one of [51], only its implementation is different as it involves a single polynomial kernel instead of a different Taylor basis and coefficients in each quadrature point on $\bar{\Gamma}_D$. Stability and error estimates can be found in [51].

*For the Aubin method:* The use of the polynomial correction in the embedded Dirichlet terms leads to

$$ 
\begin{aligned}
a\langle u, v\rangle_{\bar{\Gamma}_D} &= \gamma^{-1}(u(x), v)_{\bar{\Gamma}_D} - (\nabla u \cdot \bar{n}, v)_{\bar{\Gamma}_D}, \\
b\langle v\rangle_{\bar{\Gamma}_D} &= \gamma^{-1}(u_D, v)_{\bar{\Gamma}_D}.
\end{aligned} 
\quad (42) $$

### 5.1.2 Polynomial correction on the Neumann boundary

Neumann boundaries are treated in a similar way. However, in this case, the boundary condition also depends on the normal, which must be accounted for in the correction. We consider the vectors

$$ \bar{n} = \bar{n}(\bar{x}), \quad n = n(x) = n(x), \quad t = t(x) = t(x), \quad (43) $$

that represents the unit normal vector to the surrogate boundary and the corresponding perpendicular unit normal and tangent vectors to the true boundary. We then introduce the decomposition $\bar{n} = n(\bar{n} \cdot n) + t(\bar{n} \cdot t)$ and write that

$$ \bar{q}_N = (\bar{n} \cdot n)\nabla u(\bar{x}) \cdot n + (\bar{n} \cdot t)\nabla u(\bar{x}) \cdot t. \quad (44) $$

In the second term, we modify using the true Neumann data. As for the Dirichlet conditions, e.g., see [39], we introduce a Taylor series of the gradient in the fixed normal direction, $n$, which is matched with the true Neumann data as

$$ q_N(x) = \nabla u(\bar{x}) \cdot n + D_N^k, \quad \text{where} \quad D_N^k = \sum_{\alpha=1}^k \frac{1}{\alpha!} (\nabla^\alpha) (\nabla u \cdot n)(\bar{x}) d(\bar{x})^\alpha. \quad (45) $$

As discussed in [50], one should be careful about how the derivative terms in the expansion above are obtained. If they are evaluated from a finite element cell polynomial of degree $k$, the $k$'th derivative is not consistent. Unless something special is done to recover this, the above development can only be performed to degree $k-1$. We use the identity

$$ D_N^k = \nabla u(x) \cdot n - \nabla u(\bar{x}) \cdot n, \quad (46) $$

and combine it with (44). After a few manipulations, we obtain

$$ \bar{q}_N = \nabla u(\bar{x}) \cdot \bar{n} + (\bar{n} \cdot n)(q_N(x) - \nabla u(x) \cdot n). \quad (47) $$

Replacing this simple expression in the bilinear form and rearranging, we get

$$ 
\begin{aligned}
a\langle u, v\rangle_{\bar{\Gamma}_N} &= -(\nabla u(\bar{x}) \cdot \bar{n}, v)_{\bar{\Gamma}_N} + (\nabla u(x) \cdot n, (\bar{n} \cdot n)v)_{\bar{\Gamma}_N}, \\
b\langle v\rangle_{\bar{\Gamma}_N} &= (q_N, (\bar{n} \cdot n)v)_{\bar{\Gamma}_N}.
\end{aligned} 
\quad (48) $$

This modification corresponds to a correction term based on the difference between the Neumann data and the cell normal gradient of the finite element polynomial evaluated on the true boundary $\Gamma$. Note that in the current formulation, we obtain a sub-optimal method with consistency at most $O(h^k)$ in the $L^2$ norm for an approximation of degree $k$. This is linked to the approximation accuracy of the first derivatives of the finite element polynomial.

## 6 Asymptotic preserving Robin formulations: one implementation for all embedded boundaries

For the Robin conditions, we examine the idea of constructing AP approximations that can be used for a single high-order implementation for all types of boundary conditions. Thus, we focus on the full Robin problem in (2), and start from the conformal weak formulation in Nitsche form given by (19). With the notation of Section 4.5, we define the Dirichlet form and Neumann form (with added penalty as in Proposition 4.2) as

$$ D := \mathcal{I}_{\Omega}(u) - (\partial_n u, v)_{\Gamma} - (f, v)_{\Omega} + \mathcal{N}_{\Gamma}(u - u_{RD}), \quad (49) $$
$$N := \mathcal{I}_{\Omega}(u) - (\partial_n u, v)_{\Gamma} - (f, v)_{\Omega} + \mathcal{N}_{\Gamma}(\partial_n u - q_{RN}). \quad (50)$$

As discussed in [57] it is easy to see that (19) can be obtained as the linear combination

$$\frac{\gamma}{\gamma + \varepsilon} D + \frac{\varepsilon}{\gamma + \varepsilon} N = 0. \quad (51)$$

It is quite tempting to perform a similar combination for the embedded case to obtain a SBM Robin formulation. Unfortunately, this may lead to an inconsistent method. To see this, let us set

$$\mathcal{N}_{\Gamma}(u) := (u, v(x)/\gamma - \partial_{\bar{n}} v(\bar{x})). \quad (52)$$

Using (41) and (48), we define the embedded Robin problem from

$$\bar{D} := \mathcal{I}_{\bar{\Omega}}(u) - (\partial_{\bar{n}} u(\bar{x}), v)_{\bar{\Gamma}} - (f, v)_{\bar{\Omega}} + \mathcal{N}_{\bar{\Gamma}}(u(\bar{x}) - u_{RD}), \quad (53)$$

and

$$\bar{N} := \mathcal{I}_{\bar{\Omega}}(u) - (\partial_{\bar{n}} u(\bar{x}), v)_{\bar{\Gamma}} - (f, v)_{\bar{\Omega}} + \gamma \mathcal{N}_{\bar{\Gamma}}((\bar{n} \cdot n)(\partial_n u(\bar{x}) - q_{RN})). \quad (54)$$

Quick computations show that the method defined by

$$\frac{\gamma}{\gamma + \varepsilon} \bar{D} + \frac{\varepsilon}{\gamma + \varepsilon} \bar{N} = 0, \quad (55)$$

leads to the variational formulation

$$\mathcal{I}_{\bar{\Omega}}(u) - (\partial_{\bar{n}} u(\bar{x}), v)_{\bar{\Gamma}} - (f, v)_{\bar{\Omega}} + \frac{\gamma}{\gamma + \varepsilon} \mathcal{N}_{\bar{\Gamma}}(u(\bar{x}) - u_{RD}) + \frac{\gamma \varepsilon}{\gamma + \varepsilon} \mathcal{N}_{\bar{\Gamma}}((\bar{n} \cdot n)(\partial_n u(\bar{x}) - q_{RN})) = 0. \quad (56)$$

By construction, for $\varepsilon = 0$ this method reduces exactly to (41), and – after some manipulations – for $\varepsilon^{-1} \to 0$ to

$$\mathcal{I}_{\bar{\Omega}}(u) - (\partial_{\bar{n}} u(\bar{x}), v)_{\bar{\Gamma}} - (f, v)_{\bar{\Omega}} + (\partial_n u(\bar{x}) - q_{RN}, (\bar{n} \cdot n)v)_{\bar{\Gamma}} - \gamma(\partial_n u(\bar{x}) - q_{RN}, (\bar{n} \cdot n)\partial_{\bar{n}}v)_{\bar{\Gamma}} = 0, \quad (57)$$

which is equal to (48) with the extra penalty also found in (22). However, if we compute the Euler-Lagrange equations for general values of $\varepsilon$, we obtain

$$(-\nabla^2 u + \alpha u - f, v)_{\bar{\Omega}} + \frac{\gamma}{\varepsilon + \gamma} \underbrace{(u(\bar{x}) - u_{RD} + \varepsilon(\bar{n} \cdot n)(\partial_n u(\bar{x}) - q_{RN}))}_{\bar{N}} v(\bar{x})/\gamma - \partial_{\bar{n}} v(\bar{x})_{\bar{\Gamma}} = 0, \quad (58)$$

which in general is not consistent, even within a given power of $h$, due to the presence of the $(\bar{n} \cdot n)$ product in the under-braced term that depends strongly on the generated mesh topology and true domain shape.

In the following sections, we propose three approaches that correct this problem. These methods will provide Robin formulations that are not only consistent but also enjoy an AP property. Hereby, they can be used to approximate any boundary condition simply by modifying the value of $\varepsilon$.

## 6.1 The Nitsche-Robin method based on corrected linear combination coefficients

We first consider a simple correction of the coefficients in the linear combination discussed in the previous section. As with the inconsistent approach, we use (39) and (47) to define $\bar{u}_{RD}$ and $\bar{q}_{RN}$, given $u_{RD}$ and $q_{RN}$. However, we now set $\bar{\gamma} = \bar{n}_n \gamma$, by introducing the short notation

$$\bar{n}_n := \bar{n} \cdot n. \quad (59)$$

We now consider the modified variant of (56) that reads

$$\mathcal{I}_{\bar{\Omega}}(u) - (\partial_{\bar{n}} u(\bar{x}), v)_{\bar{\Gamma}} - (f, v)_{\bar{\Omega}} + \mathcal{N}_{\bar{\Gamma}} \left( \frac{\bar{\gamma}}{\bar{\gamma} + \varepsilon} (u(\bar{x}) - u_{RD}) \right) + \mathcal{N}_{\bar{\Gamma}} \left( \frac{\bar{\gamma} \varepsilon}{\bar{\gamma} + \varepsilon} (\partial_n u(\bar{x}) - q_{RN}) \right) = 0. \quad (60)$$

The method is equivalent to the combination

$$\frac{\bar{\gamma}}{\bar{\gamma} + \varepsilon} \bar{D} + \frac{\varepsilon}{\bar{\gamma} + \varepsilon} \bar{N} = 0, \quad (61)$$

which, however, due to the variability of $\bar{n}_n$, only makes sense when the combination coefficients are finally included in the boundary integrals as done in (60). A quick check shows that the limits for $\varepsilon \to 0$ and $\varepsilon^{-1} \to 0$ are the same as (56). Moreover, the Euler-Lagrange equations in the general case read

$$(-\nabla^2 u + \alpha u - f, v)_{\bar{\Omega}} + \underbrace{(u(\bar{x}) - u_D + \varepsilon(\partial_n u(\bar{x}) - q_{RN}))}_{\bar{N}} \frac{\bar{\gamma}}{\bar{\gamma} + \varepsilon} (v(\bar{x})/\gamma - \partial_{\bar{n}} v(\bar{x}))_{\bar{\Gamma}} = 0, \quad (62)$$

where the under-braced term is consistent with the full Robin boundary condition within the accuracy of the extrapolation. For smooth enough problems, the latter is expected to behave as $\max(\mathcal{O}(h^{k+1}), \min(\varepsilon, C)\mathcal{O}(h^k))$ for some $C > 0$. More interestingly, the above approximation is asymptotic preserving, as shown below.
**Proposition 6.1** (AP property of the corrected coefficients Nitsche-Robin embedded method). *The Robin embedded method with shifted boundary polynomial correction defined by (60) is asymptotic preserving in both the Dirichlet limit $\epsilon \to 0$ and in the Neumann limit $\epsilon^{-1} \to 0$. In particular, we have that*

* in the limit $\epsilon \ll 1$, the discrete modes {$u_m$}$_{m \ge 0}$ expanded as in (3a) verify the coercive asymptotic variational problems

$$ \mathcal{I}_{\overline{\Omega}}(u_m) - (\partial_{\bar{n}} u_m(\bar{x}), v)_{\overline{\Gamma}} - (f_m, v)_{\overline{\Omega}} + \mathcal{N}_{\overline{\Gamma}}(u_m(x) - u_{D,m}) = \overline{\mathcal{E}}_m, \quad (63) $$

with $u_{D,m}$ defined by Proposition 2.1, and with $f_0 = f$ and $f_m = 0$ for $m \ge 1$, and where $\overline{\mathcal{E}}_0 = 0$ while

$$ \overline{\mathcal{E}}_m := - \sum_{\ell=0}^{m-1} (-\bar{\gamma})^{-(m-\ell)} \mathcal{N}_{\overline{\Gamma}}(u_{\ell}(\mathbf{x}) - u_{D,\ell}), \quad \forall m \ge 1. \quad (64) $$

* in the limit $\epsilon^{-1} \ll 1$, discrete modes expanded as in (4a) verify the asymptotic variational problems

$$ \mathcal{I}_{\overline{\Omega}}(u_m) - (\partial_{\bar{n}} u_m(\bar{x}), v)_{\overline{\Gamma}} - (f_m, v)_{\overline{\Omega}} + \mathcal{N}_{\overline{\Gamma}}(\bar{\gamma}(\nabla u_m(\mathbf{x}) \cdot \mathbf{n} - q_{N,m})) = \overline{\mathcal{E}}_m, \quad (65) $$

with $q_{N,m}$ defined by Proposition 2.2, and with $f_0 = f$ and $f_m = 0$ for $m \ge 1$, and where $\overline{\mathcal{E}}_0 = 0$ while

$$ \overline{\mathcal{E}}_m := \sum_{\ell=0}^{m-1} (-\bar{\gamma})^{(m+1-\ell)} \mathcal{N}_{\overline{\Gamma}}(\nabla u_{\ell}(\mathbf{x}) \cdot \mathbf{n} - q_{N,\ell}), \quad \forall m \ge 1. \quad (66) $$

*Proof.* To obtain the proof, we proceed by inserting the asymptotic developments formally into (60). To account for the effect of the $1/(\bar{\gamma} + \epsilon)$ term, we cannot multiply the whole equation by it due to its variability in space. Instead, we replace this term by using a Maclaurin series expansion as

$$ \frac{1}{1+x} = \sum_{n \ge 0} (-x)^n, \quad (67) $$

applied to $x = \epsilon\bar{\gamma}^{-1}$ in the limit $\epsilon \ll 1$, and to $\epsilon^{-1}\bar{\gamma}$ in the limit $\epsilon^{-1} \ll 1$. Both (63) and (65) are obtained classically by gathering the terms multiplying equal powers of the small parameter. Each asymptotic problem is a Nitsche-type approximation of the appropriate order continuous problem, with error terms proportional to the boundary errors. In particular, the variational forms (63) are equivalent to the symmetrized Nitsche SBM approximation of the corresponding continuous problems on non-conformal meshes, whose coercivity has been proven in [51]. □

As for the inconsistent approach in (56), this method allows us to use the surrogate data independently in the Neumann and Dirichlet cases. Yet, the method restores the consistency with the Robin problem for general values of $\epsilon$ and is asymptotic preserving. However, it involves coefficients depending on the scalar product in (59), whose value cannot be a priori set. This modifies the weighting of the boundary terms in an uncontrolled – albeit possibly small – manner.

Concerning the parameter $\gamma$, the choice $\gamma = C_\gamma h$ allows us to ensure the coercivity of the bilinear form in the Neumann asymptotic limit for the conformal case. Unfortunately, stability results are only available for pure Dirichlet boundaries, and the impact of the penalty on this property seems to be relatively important [52, 76]. Thus, the same definition of $\gamma$ – as the one given in the conformal case – seems to be the most appropriate.

For completeness, we report the bi-linear and linear forms corresponding to the corrected scheme, which read

$$ a\langle u, v \rangle_{\overline{\Gamma}_R} = -(\partial_{\bar{n}} u(\bar{x}), v)_{\overline{\Gamma}_R} + \mathcal{N}_{\overline{\Gamma}_R} \left( \frac{\bar{\gamma}}{\bar{\gamma} + \epsilon} u(\bar{x}) \right) + \mathcal{N}_{\overline{\Gamma}_R} \left( \frac{\bar{\gamma}\epsilon}{\bar{\gamma} + \epsilon} \partial_n u(\bar{x}) \right), \quad (68) $$

$$ b\langle v \rangle_{\overline{\Gamma}_R} = \mathcal{N}_{\overline{\Gamma}_R} \left( \frac{\bar{\gamma}}{\bar{\gamma} + \epsilon} u_{RD} \right) + \mathcal{N}_{\overline{\Gamma}_R} \left( \frac{\bar{\gamma}\epsilon}{\bar{\gamma} + \epsilon} q_{RN} \right). $$

## 6.2 The Nitsche-Robin method using polynomial corrections in the full Robin condition

We now proceed differently by working directly with the full Robin condition. We modify the following SBM approach discussed in Sections 5.1.1 and 5.1.2 as

$$ \begin{aligned} u - u_{RD} + \epsilon(\partial_n u - q_{RN}) &= u(\bar{x}) + \mathcal{D}_D^k - u_{RD} + \mathcal{O}(\|\mathbf{d}^{k+1}\|) + \epsilon(\partial_n u(\bar{x}) + \mathcal{D}_N^k - q_{RN}) + \epsilon\mathcal{O}(\|\mathbf{d}^k\|) \\ &= u(\bar{x}) + \mathcal{D}_D^k - u_{RD} + \mathcal{O}(\|\mathbf{d}^{k+1}\|) + \epsilon(\partial_{\bar{n}} u(\bar{x}) + \mathcal{D}_N^k + \mathcal{D}_N^n - q_{RN}) + \epsilon\mathcal{O}(\|\mathbf{d}^k\|), \end{aligned} \quad (69) $$
having denoted by $\mathcal{D}_{D/N}^k$ the high-order polynomial corrections and set

$$ \mathcal{D}_N^n := \partial_n u(\bar{x}) - \partial_{\bar{n}} u(\bar{x}). \qquad (70) $$

We now define

$$ \bar{u}_{RD} := u_{RD} - \mathcal{D}_D^k, \quad \text{and} \quad \bar{q}_{RN} := q_{RN} - \mathcal{D}_N^k - \mathcal{D}_N^n, \qquad (71) $$

and apply the generalized Nitsche's method on the surrogate domain with this data. Using the polynomial corrections based on the cell finite element expansion, we obtain the problem

$$ \mathcal{I}_{\bar{\Omega}}(u) - (\partial_{\bar{n}} u(\bar{x}), v)_{\bar{\Gamma}} - (f, v)_{\bar{\Omega}} + \frac{\gamma}{\gamma + \varepsilon} \mathcal{N}_{\bar{\Gamma}}(u(\bar{x}) - u_{RD}) + \frac{\gamma \varepsilon}{\gamma + \varepsilon} \mathcal{N}_{\bar{\Gamma}}(\partial_n u(\bar{x}) - q_{RN}) = 0. \qquad (72) $$

A few computations show that the Euler-Lagrange equations for this problem are

$$ (-\nabla^2 u + \alpha u - f, v)_{\bar{\Omega}} + \frac{\gamma}{\gamma + \varepsilon} \underbrace{(u(\bar{x}) - u_D + \varepsilon(\partial_n u(\bar{x}) - q_{RN}))}_{\bar{\Gamma}} v(\bar{x}) / \gamma - \partial_{\bar{n}} v(\bar{x})_{\bar{\Gamma}} = 0, \qquad (73) $$

which readily shows the consistency with the Robin condition within an error of order $\max(\mathcal{O}(h^{k+1}), \min(\varepsilon, C)\mathcal{O}(h^k))$ for some constant $C > 0$. The above approximation is also asymptotic preserving, and we can prove the following.

**Proposition 6.2 (AP property of the full Nitsche-Robin embedded method).** The Robin embedded method with shifted boundary polynomial correction (72) is asymptotic preserving in both the Dirichlet limit $\varepsilon \to 0$ and in the Neumann limit $\varepsilon^{-1} \to 0$. In particular, we have that

* in the limit $\varepsilon \ll 1$, the discrete modes $\{u_m\}_{m \ge 0}$ expanded as in (3a) verify the coercive asymptotic variational problems

$$ \mathcal{I}_{\bar{\Omega}}(u_m) - (\partial_{\bar{n}} u_m(\bar{x}), v)_{\bar{\Gamma}} - (f_m, v)_{\bar{\Omega}} + \mathcal{N}_{\bar{\Gamma}}(u_m(\bar{x}) - u_{D,m}) = \bar{\mathcal{E}}_m, \qquad (74) $$

with $u_{D,m}$ defined by Proposition 2.1, and with $f_0 = f$ and $f_m = 0$ for $m \ge 1$, and where $\bar{\mathcal{E}}_0 = 0$ while

$$ \bar{\mathcal{E}}_m := - \sum_{\ell=0}^{m-1} (-\gamma)^{-(m-\ell)} \mathcal{N}_{\bar{\Gamma}}(u_{\ell}(\bar{x}) - u_{D,\ell}), \quad \forall m \ge 1. \qquad (75) $$

* in the limit $\varepsilon^{-1} \ll 1$, discrete modes expanded as in (4a) verify the asymptotic variational problems

$$ \mathcal{I}_{\bar{\Omega}}(u_m) - (\partial_{\bar{n}} u_m(\bar{x}), v)_{\bar{\Gamma}} - (f_m, v)_{\bar{\Omega}} + \mathcal{N}_{\bar{\Gamma}}(\gamma(\nabla u_m(\bar{x}) \cdot \mathbf{n} - q_{N,m})) = \bar{\mathcal{E}}_m, \qquad (76) $$

with $q_{N,m}$ defined by Proposition 2.2, and with $f_0 = f$ and $f_m = 0$ for $m \ge 1$, and where $\bar{\mathcal{E}}_0 = 0$ while

$$ \bar{\mathcal{E}}_m := \sum_{\ell=0}^{m-1} (-\gamma)^{(m+1-\ell)} \mathcal{N}_{\bar{\Gamma}}(\nabla u_{\ell}(\bar{x}) \cdot \mathbf{n} - q_{N,\ell}) \quad \forall m \ge 1. \qquad (77) $$

*Proof.* The proof follows the steps of the proof of Proposition 4.1 and 4.2. It is omitted for brevity. □

This approach has the advantage of not involving a parameter depending on the orientation of the surrogate boundary compared to the true boundary. For completeness, we report the bi-linear and linear forms corresponding to the corrected scheme, which read

$$ a\langle u, v \rangle_{\bar{\Gamma}_R} = - (\partial_{\bar{n}} u(\bar{x}), v)_{\bar{\Gamma}_R} + \frac{\gamma}{\gamma + \varepsilon} \mathcal{N}_{\bar{\Gamma}_R}(u(\bar{x})) + \frac{\gamma \varepsilon}{\gamma + \varepsilon} \mathcal{N}_{\bar{\Gamma}_R}(\partial_n u(\bar{x})), \qquad (78) $$

$$ b\langle v \rangle_{\bar{\Gamma}_R} = \frac{\gamma}{\gamma + \varepsilon} \mathcal{N}_{\bar{\Gamma}_R}(u_{RD}) + \frac{\gamma \varepsilon}{\gamma + \varepsilon} \mathcal{N}_{\bar{\Gamma}_R}(q_{RN}). \qquad (79) $$

## 6.3 The Aubin-Robin formulation

We consider the embedded formulation of the Aubin-Robin obtained by simplifying the last approach. In particular, using the notation

$$ \mathcal{A}_{\bar{\Gamma}}(u) := (u, v(\bar{x}) / \gamma)_{\bar{\Gamma}}, \qquad (79) $$

we consider the embedded Aubin-Robin method obtained setting as

$$ a\langle u, v \rangle_{\bar{\Gamma}_R} = - (\partial_{\bar{n}} u(\bar{x}), v)_{\bar{\Gamma}_R} + \frac{\gamma}{\gamma + \varepsilon} \mathcal{A}_{\bar{\Gamma}_R}(u(\bar{x})) + \frac{\gamma \varepsilon}{\gamma + \varepsilon} \mathcal{A}_{\bar{\Gamma}_R}(\partial_n u(\bar{x})), \qquad (80) $$

$$ b\langle v \rangle_{\bar{\Gamma}_R} = \frac{\gamma}{\gamma + \varepsilon} \mathcal{A}_{\bar{\Gamma}_R}(u_{RD}) + \frac{\gamma \varepsilon}{\gamma + \varepsilon} \mathcal{A}_{\bar{\Gamma}_R}(q_{RN}). \qquad (81) $$

This method can be shown to be AP, and in particular, enjoy the following property.
**Proposition 6.3** (AP property of the Aubin-Robin embedded method). The Robin embedded method with shifted boundary polynomial correction (72) is asymptotic preserving in both the Dirichlet limit $\epsilon \to 0$ and in the Neumann limit $\epsilon^{-1} \to 0$. In particular, we have that

* in the limit $\epsilon \ll 1$, the discrete modes $\{u_m\}_{m \ge 0}$ expanded as in (3a) verify the coercive asymptotic variational problems

$$ \mathcal{I}_{\bar{\Omega}}(u_m) - (\partial_{\bar{n}} u_m(\bar{x}), v)_{\bar{\Gamma}} - (f_m, v)_{\bar{\Omega}} + \mathcal{A}_{\bar{\Gamma}}(u_m(\bar{x}) - u_{D,m}) = \bar{\mathcal{E}}_m, \quad (81) $$

with $u_{D,m}$ defined by Proposition 2.1, and with $f_0 = f$ and $f_m = 0$ for $m \ge 1$, and where $\bar{\mathcal{E}}_0 = 0$ while

$$ \bar{\mathcal{E}}_m := - \sum_{\ell=0}^{m-1} (-\gamma)^{-(m-\ell)} \mathcal{A}_{\bar{\Gamma}}(u_{\ell}(\mathbf{x}) - u_{D,\ell}), \quad \forall m \ge 1. \quad (82) $$

* in the limit $\epsilon^{-1} \ll 1$, discrete modes expanded as in (4a) verify the asymptotic variational problems

$$ \mathcal{I}_{\bar{\Omega}}(u_m) - (\partial_{\bar{n}} u_m(\bar{x}), v)_{\bar{\Gamma}} - (f_m, v)_{\bar{\Omega}} + \mathcal{A}_{\bar{\Gamma}}(\gamma(\nabla u_m(\bar{x}) \cdot \mathbf{n} - q_{N,m})) = \bar{\mathcal{E}}_m, \quad (83) $$

with $q_{N,m}$ defined by Proposition 2.2, and with $f_0 = f$ and $f_m = 0$ for $m \ge 1$, and where $\bar{\mathcal{E}}_0 = 0$ while

$$ \bar{\mathcal{E}}_m := \sum_{\ell=0}^{m-1} (-\gamma)^{(m+1-\ell)} \mathcal{A}_{\bar{\Gamma}}(\nabla u_{\ell}(\mathbf{x}) \cdot \mathbf{n} - q_{N,\ell}) \quad \forall m \ge 1. \quad (84) $$

*Proof.* The proof follows the steps of the proof of Proposition 4.1 and 4.2. It is omitted for brevity. □

## 6.4 Choice of the embedding: a note on interpolation error and conditioning

It is well known that the condition number of the system matrix, $\kappa(\mathcal{A})$, is of immense importance in relation to the cost of the matrix inversion and the accuracy of the resulting solution. Important factors for the value of $\kappa(\mathcal{A})$ are the choice of the finite element basis and associated quadrature points. The iso-parametric SEM is highly optimized for well-behaved interpolation, differentiation, and integration within the element at the delicately picked quadrature points [18, 71]. However, extrapolation may be as important as interpolation when considering unfitted finite elements. In the SBM setting using the classical extrapolation of $u(\mathbf{x})$ and $\nabla u(\mathbf{x})$ will therefore result in evaluating the finite element polynomial outside the local element.

It is also well known that the Lebesque constant, $\Lambda_P$, is a good measure of the polynomial interpolation quality, e.g., see [77]. The Lebesque constant is defined by the Lebesque function, $\mathcal{L}_P = \mathcal{L}_P(\mathbf{r})$, in the reference domain, $\mathcal{R}$, as

$$ \Lambda_P = \max_{\mathbf{r} \in \mathcal{R}} \mathcal{L}_P(\mathbf{r}), \quad \text{where} \quad \mathcal{L}_P(\mathbf{r}) = \sum_{i=1}^{N_{ep}} |h_i(\mathbf{r})|. \quad (85) $$

From [78], we see that the interpolation error in the maximum norm, $\|f - \mathcal{I}_P f\|_{\infty}$, is bounded from above by the Lebesque constant as

$$ \|f - \mathcal{I}_P f\|_{\infty} \le (1 + \Lambda_P) \|f - f^*\|_{\infty}, \quad (86) $$

<figure><img src="image_5.png" /><figcaption>Figure 3: Graphical setup for the numerical domain when computing the Lebesque constants in interpolation on $\mathcal{R}$ and extrapolation on $C - \mathcal{R}$.</figcaption></figure>
<table><caption>Table 1: Lebesque constant, $\Lambda_P$, for polynomial orders, $P = \{1, ..., 10\}$, in interpolation and extrapolation as defined on Figure 3 using the node distribution due to [71].</caption><thead><tr><th>P</th><th>Interpolation</th><th>Extrapolation</th></tr></thead><tbody><tr><td>1</td><td>1.00 · 10<sup>0</sup></td><td>2.76 · 10<sup>0</sup></td></tr><tr><td>2</td><td>1.67 · 10<sup>0</sup></td><td>1.42 · 10<sup>1</sup></td></tr><tr><td>3</td><td>2.11 · 10<sup>0</sup></td><td>8.03 · 10<sup>1</sup></td></tr><tr><td>4</td><td>2.66 · 10<sup>0</sup></td><td>4.67 · 10<sup>2</sup></td></tr><tr><td>5</td><td>3.12 · 10<sup>0</sup></td><td>2.73 · 10<sup>3</sup></td></tr><tr><td>6</td><td>3.70 · 10<sup>0</sup></td><td>1.60 · 10<sup>4</sup></td></tr><tr><td>7</td><td>4.27 · 10<sup>0</sup></td><td>9.44 · 10<sup>4</sup></td></tr><tr><td>8</td><td>4.96 · 10<sup>0</sup></td><td>5.59 · 10<sup>5</sup></td></tr><tr><td>9</td><td>5.74 · 10<sup>0</sup></td><td>3.32 · 10<sup>6</sup></td></tr><tr><td>10</td><td>6.67 · 10<sup>0</sup></td><td>1.99 · 10<sup>7</sup></td></tr></tbody></table>

where $f$ is some function with the interpolation $\mathcal{I}_P f$ over $\mathcal{R}$ and $f^* \in \mathbb{P}^P$. The nodal Lagrange basis is bounded within the reference domain, yet, the basis functions grow unbounded outside at an exponential rate as higher and higher orders are considered. Hereby, it is expected that extrapolating the basis function will increase the Lebesque constant such that the interpolation error increases.

To support this argument, a simple numerical exercise is considered inspired by [77]. The results are highlighted in Table 1, where the Lebesque constants are computed numerically for the classical interpolation case on $\mathcal{R}$ and similarly computed in extrapolation. The latter case is done for a circle of radius $R = 1.75$ with its center at the barycenter of the reference triangle, $(r, s) = -1/3$, i.e., $C = \{r = (r, s)|(r + 1/3)^2 + (s + 1/3)^2 \le 1.75^2\}$. The extrapolation is done on $C - \mathcal{R}$. This setup is visualized in Figure 3. From the Table, it is clear that the Lebesque constant – and therefore also the extrapolation error – grows rapidly and unbounded when increasing the polynomial orders. Note that the tabulated values for interpolation agree with the ones presented in [70, 77].

This motivates a comparison of the extrapolation-based SBM thoroughly with the one using the interpolation domain, as introduced in Section 3.1. In particular, in the upcoming, we will compare the following four approaches denoted by the abbreviations: **CBM**, **SBM-e**, **SBM-ei**, and **SBM-i** as

1. **CBM**: Body-fitted affine iso-parametric meshes where the elements conform to $\Gamma$. See Figure 4 (top-left).

2. **SBM-e**: Unfitted affine meshes, where $\bar{\Omega}$ is enclosed completely within $\Gamma$, and the mapping is defined in the original SBM fashion, i.e., along the true normal, $n$, on $\Gamma$. This approach results in a pure extrapolation procedure as seen in Figure 4 (top-right). Note how all points are in pure extrapolation.

3. **SBM-ei**: Unfitted meshes, where $\Gamma$ is within $\bar{\Omega}$, and the mapping is along $n$ on $\Gamma$ as for the SBM-e. As visible in Figure 4 (lower-left), not all the quadrature points are mapped into the cells, thus resulting in a mixed extrapolation and interpolation procedure.

4. **SBM-i**: Unfitted meshes, where $\Gamma$ is within $\bar{\Omega}$, and the mapping is defined locally such that any point $\bar{x}$ of a surrogate boundary edge is mapped to a point $x = M(\bar{x})$ within the element to which the edge belongs, as depicted in Figure 4 (lower-right). For this paper, we distribute the points on $\Gamma$ equidistantly between the edge intersections.

# 7 Numerical Experiments

In the following, we show a number of numerical experiments with Dirichlet, Neumann, and Robin boundary conditions. These will be of the Aubin and Nitsche variational forms using the four boundary and mapping approaches (CBM, SBM-e, SBM-ei, and SBM-i). See Figure 4 for the general concept hereof. We consider meshes composed of affine triangular elements that are quasi-uniform in size and arranged in an unstructured manner. Finite elements up to polynomial order $P = 10$ are used on meshes that are generated with the open-source meshing software Gmsh [79]. An essential input to Gmsh is the characteristic length, $l_c$, which is an approximate size of the desired elements, $l_c \approx h_{avg}$. The true geometrical domain, $\Omega$, is of a circle with radius $R = 0.375$, unless stated otherwise. The circle is placed with center at $(x_0, y_0) = (L_x/2, L_y/2)$, with $L_x = 1$ and $L_y = 1$ (again, unless stated otherwise) being the domain lengths of the background mesh in the $x$ and $y$-directions, respectively. The proposed model is currently implemented in Matlab, version R2022a, and the linear system of equations in (17) is solved with a direct solver within this framework.
<figure><img src="image_6.png" /><figcaption>Figure 4: Illustrations of CBM (top-left) and the three SBM procedures, SBM-e (top-right), SBM-ei (lower-left), and SBM-i (lower-right), for $P = 3$. The red line is $\Gamma$, the blue lines constitute the elements that form $\bar{\Omega}_h$, the red and blue dots represent nodal points $x \in \Gamma$ and $\bar{x} \in \bar{\Gamma}$, respectively, and the dashed black line is the mapping, $M(\bar{x})$.</figcaption></figure>

For simplicity, we use the same exact solution throughout the paper obtained by the classical method of manufactured solutions (MMS). We construct the solution to the Poisson problem, such that $u_{\text{MMS}}(x, y) \in C^{\infty}$, by combining trigonometric functions with a linear polynomial twist as

$$u_{\text{MMS}}(x, y) = \cos\left(5\frac{x\pi}{L_x}\right) \sin\left(5\frac{y\pi}{L_y}\right) + 2x - 1y. \qquad (87)$$

As an error measure, we exploit exact global integration through the SEM on the affine elements, with this $\|u_h - u\|_{L^1} = \int_{\Omega} |u_{\text{SEM}} - u_{\text{MMS}}|d\Omega$. Also, we opt to use a second measure as $\|Au - b\|_{L^1}$ where $u = u_{\text{MMS}}$. This measure returns the level of truncation error without the influence of, e.g., conditioning issues from the matrix inversion.

**Outline of experiments:** We consider the three types of boundary conditions one by one. For Dirichlet, a verification of Aubin's and Nitsche's formulations is carried out on aligned meshes. This numerical experiment aims to – qualitatively – compare with other high-order SBM implementations in the literature. Then, we consider the solution on completely unstructured meshes, where $(\bar{n} \cdot n)$ is not necessarily close to unity. For Neumann conditions, we limit ourselves – for brevity – to a single p-convergence study; however, general comments for the Neumann problem are included. For the Robin formulation, the consistent and inconsistent formulations are compared for Aubin and Nitsche. Also included is an assessment on unstructured meshes and high-order verification of the implementation. Lastly, we use the Robin condition with variable $\epsilon$ to solve a mixed Dirichlet-Neumann problem and confirm the asymptotic preserving properties.
<figure><img src="image_7.png" /><figcaption>Figure 5: Mesh visualization for overall assessment. 100 circles of radius $R = 0.375$ with randomly generated centers. Highlighted center in red and all centers in blue.</figcaption></figure>

## 7.1 Dirichlet boundary conditions

First, we consider the pure Dirichlet problem, where we choose the penalty constant to be $\gamma^{-1} = 2/h_{\text{avg}}$, as this – in the CBM case – has provided the expected algebraic convergence rates in general. The authors note that the choice of an insufficient penalty constant can limit the convergence, yet too large values can amplify potential issues with ill-conditioning. A penalty constant that varies with the approximation order is presented in [51], whereas a penalty-free high-order method is demonstrated in [52].

### 7.1.1 Verification on aligned meshes

A $h$-convergence comparison is carried out for $P = \{1, ..., 5\}$ for Aubin and Nitsche using CBM, SBM-e, SBM-ei, and SBM-i. In Appendix A, the figures are presented for Aubin in Figure 16 and Nitsche in Figure 17. We purposely align the surrogate boundary to the true boundary, such that $\bar{n} \cdot n \approx 1$. To do so, we generate the unfitted mesh by fitting another circle of radius $R \pm l_c/4$. The cases with "+" are for interpolation (SBM-ei and SBM-i), whereas "-" is for extrapolation (SBM-e). Meshes for CBM are constructed trivially. Alongside the $h$-convergence study, the associated matrix conditioning, $\kappa(\mathcal{A})$, is considered as well. As shown in all previous works, despite unfitted meshes, the shifted boundary approach genuinely provides high-order results with comparable errors. All results show optimal converged behavior of $O(h_{\text{max}}^{P+1})$ and align – in magnitude and scaling – with those presented in [51, 52]. Ultimately, this confirms the potential of the SBM strategy. The condition number obtained is comparable in size where all the methods somewhat follow monotonically the $O(1/h_{\text{max}}^2)$ scaling as for the CBM approach. This is one of the key features of the SBM approach due to the complete absence of small cut cells. Moreover, the Aubin formulation looks to provide slightly more regular and smooth rates compared to the Nitsche formulation.

### 7.1.2 Assessment on unstructured meshes

When considering any of the shifted boundary formulations, we find some irregularities in the convergence properties occasionally. Those irregularities are due to mesh sensitivity, i.e., where and how the true domain is embedded in the background mesh. Here, the distance between the true and surrogate boundaries can be small or larger; the boundaries can be almost parallel or almost perpendicular, etc. Due to this, the convergence rates are not always perfectly monotonic, as seen multiple times in SBM literature [50, 52]. Therefore, we wish to compare SBM-e, SBM-ei, and SBM-i more quantitatively to justify numerically whether the interpolating approaches improve the numerical scheme over the classical extrapolation approach. For this, a $(L_x \times L_y)$ square domain is considered with $L_x = L_y = 2$ which is meshed using $l_c = 0.15$. Now, 100 circles are embedded into the background mesh (one at a time) given a randomly generated center as shown in Figure 5.

For each of the 100 mesh configurations, we consider $p$-convergence studies and the associated matrix conditioning. The median and the variance of $||Au - b||_{L^1}$ and $\kappa(\mathcal{A})$ of the 100 discrete solutions are then computed and displayed
<figure><img src="image_8.png" /><figcaption>Figure 6: Median (top row) and variance (bottom row) of 100 discrete solutions with the setup as in Figure 5 for the Dirichlet problem. Quantities are for $\|Au - b\|_{L^1}$ and $\kappa(A)$ for Aubin and Nitsche formulation with SBM-e, SBM-ei, and SBM-i.</figcaption></figure>

in 6 for SBM-e, SBM-ie, and SBM-e with both Aubin and Nitsche formulations. Regarding the median, we observe for the Aubin formulation how SBM-ei and SBM-i provide a more favorable convergence rate than SBM-e. A similar conclusion is present for the Nitsche formulations; however, SBM-i looks to be slightly advantageous over SBM-ei. The same observation is seen in the conditioning median, especially for $P > 3$. The differences between interpolation and extrapolation are comparable to the tabulated values in Table 1. Also, it can be noted that as long as SBM-i is used, there is a minor difference between Aubin and Nitsche in convergence and conditioning. Regarding the variances, SBM-e generally provides the largest compared to the two interpolation approaches. This is true for both convergence and conditioning. Ultimately, the results emphasize how both the Aubin and Nitsche formulations favor using interpolation over extrapolation of the basis function. The authors note that the presented numerical analysis is based on one single mesh; thus, it is not general proof. However, we have not seen any contradictory results in the research for this paper.

## 7.2 Neumann boundary conditions

Next, the pure Neumann problem is considered. For conciseness, we only present one numerical study, namely p-convergence on an aligned mesh with $l_c = 1/16$ using SBM-e and SBM-i on two boundary formulations:

* Formulation 1: Standard form.

$$a\langle u, v \rangle_{\overline{\Gamma}_N} = 0, \quad \text{and} \quad b\langle v \rangle_{\overline{\Gamma}_N} = (\bar{q}_N, v)_{\overline{\Gamma}_N}$$

* Formulation 2: Standard form with symmetrical term.

$$a\langle u, v \rangle_{\overline{\Gamma}_N} = -\gamma^{-1}(\nabla u \cdot \mathbf{n}, \nabla v \cdot \mathbf{n})_{\overline{\Gamma}_N}, \quad \text{and} \quad b\langle v \rangle_{\overline{\Gamma}_N} = (\bar{q}_N, v)_{\overline{\Gamma}_N} - \gamma^{-1}(\bar{q}_N, \nabla v \cdot \mathbf{n})_{\overline{\Gamma}_N}.$$
<figure><img src="image_9.png" /><figcaption>Figure 7: p-convergence and associated condition number study of the Neumann problem on a mesh with $l_c = 1/16$. The analysis is of two boundary formulations, Formulation 1 and Formulation 2, for SBM-e and SBM-i.</figcaption></figure>

The reasoning behind considering these two formulations is to highlight the behavior of the upcoming Robin Aubin form (Formulation 1) and Robin Nitsche form (Formulation 2) in the Neumann limit as $\epsilon^{-1} \to 0$. The results are shown in Figure 7, where it can be observed that all combinations yield spectral convergence on the mesh. It can be noted that the convergence rate is ruined for Formulation 2 using SBM-e, which is due to ill-conditioning. Moreover, a decrease in convergence rate is seen when using SBM-e over SBM-i. Also, the condition number is increased when: i) using Formulation 2 over Formulation 1 and ii) using SBM-e over SBM-i.

For conciseness, further results for the Neumann problem have been intentionally left out. We stress that a sub-optimal $O(h_{\max}^P)$ convergence rate is obtained for all SBM-e, SBM-ei, and SBM-i approaches. For CBM, optimal rate of $O(h_{\max}^{P+1})$ is reached. This sub-optimality is expected – as already discussed in Section 5.1.2 – due to the polynomial correction. Some embedded formulations that allow to recover the optimal rates are proposed, e.g., in [37, 50, 80]. Also, a general increase in conditioning is seen compared to the Dirichlet problem. The explanation hereof is due to gradients in the basis functions as required in the Neumann formulations. The SEM basis functions are defined inside the elements and grow rapidly and unbounded outside the element; hence, the gradients will enhance this polynomial behavior even more outside the element. Ultimately, this will imply an increase in the conditioning. In general, SBM-e provides larger $\kappa(\mathcal{A})$ (multiple orders of magnitude) than SBM-ei and SBM-i.

## 7.3 Robin boundary conditions

The last problem to consider is the pure Robin formulation. For this, we set $\gamma^{-1}$ as with the Dirichlet problem.

### 7.3.1 Consistent versus inconsistent formulations

As outlined in Section 6, the off-the-selves Aubin and Nitsche formulations for embedded Robin boundary conditions are inconsistent if the standard polynomial corrections are applied. Therefore, we proposed a consistent way of formulating both variational statements. Ultimately, we are left with the following inconsistent and consistent equations on the boundary derived from the Euler-Lagrange equations as

$$ \text{Inconsistent: } u(x) - u_{RD} + \varepsilon(\bar{n} \cdot n)(\nabla u(x) \cdot n - q_{RN}) = 0. $$

$$ \text{Consistent: } u(x) - u_{RD} + \varepsilon(\nabla u(x) \cdot n - q_{RN}) = 0. $$

Recall that $x = M(\bar{x})$. Now, $u_{RD}$ and $q_{RN}$ are constructed from the analytical MMS solution. Therefore, we will not be able to locate the consistency issue as the error is reduced for $(\nabla u(x) \cdot n - q_{RN})$. To circumvent this, we perturb the Robin boundary condition for (2) with an arbitrary scalar value, $\delta$, as

$$ u + \varepsilon(\nabla u \cdot n) = u_{RD} - \delta + \varepsilon(q_{RN} + \delta/\varepsilon), \quad \text{on } \Gamma_R. \qquad (88) $$

With this, the consistent and inconsistent boundary equations give

$$ \text{Inconsistent: } u(x) - u_{RD} + \varepsilon(\bar{n} \cdot n)(\nabla u(x) \cdot n - q_{RN}) = \delta(\bar{n} \cdot n - 1). $$

$$ \text{Consistent: } u(x) - u_{RD} + \varepsilon(\nabla u(x) \cdot n - q_{RN}) = \delta - \delta = 0. $$
<figure><img src="image_10.png" /><figcaption>Figure 8: Consistent formulation (full lines) versus inconsistent formulation (dotted-dashed lines). For both Aubin (left) and Nitsche (right) on four different meshes, $l_c = \{1/4, 1/8, 1/16, 1/32\}$, with $p$-convergence on each. Results using SBM-i on aligned meshes with $\varepsilon = \delta = 1$.</figcaption></figure>

From this, we see that the inconsistency error – in a weak sense – becomes $\delta(\bar{n} \cdot n - 1)$, and that the consistent formulation still satisfies the $\delta$ perturbation as expected.

In Figure 8, a $p$-convergence study is carried for $P = \{1, \dots, 10\}$ on four different meshes, $l_c = \{1/4, 1/8, 1/16, 1/32\}$. Here, both the consistent and inconsistent formulations for Aubin and Nitsche are used with $\varepsilon = \delta = 1$. Ultimately, the results show how the truncation error stalls at certain levels for both of the inconsistent formulations. These levels are of first order in relation to the mesh size, i.e., related to $(\bar{n} \cdot n)$ that approach unity for $l_c \to 0$, as expected. Both of the consistent formulations can be seen to converge properly with exponential rate to machine precision. Now that the inconsistency is located we opt to solely use the consistent formulations for the remaining part of the paper.

### 7.3.2 Assessment and verification on unstructured meshes

We performed the same statistical assessment as for the pure Dirichlet problem (recall Figure 5). For this, we use $\varepsilon = 1$, and the consistent Robin formulations. The analysis can be seen in Figure 18 in Appendix B, ultimately confirming the observations and conclusions from the Dirichlet case.

Next, a verification of the Robin implementation – in terms of $h$-convergence studies – is to be carried out. For this, we select three different values of $\varepsilon = \{10^{-1}, 10^{0}, 10^{1}\}$ and run the analysis for $P = \{2, 6, 10\}$ for both Aubin and Nitsche. We compare the SBM-i approach to the reference CBM solution. Recall that the latter converges optimally in both the Dirichlet and Neumann limits, whereas the polynomial corrected SBM methods converge sub-optimal in the Neumann limit and optimal in the Dirichlet limit. The results can be seen in Figure 9. From the figure, it can be observed how the SBM-i genuinely agrees with the CBM rates, especially in the low $\varepsilon$ case. Moreover, as $\varepsilon$ increases, the difference in convergence rates (optimal versus sub-optimal) appears. Lastly, the authors note the slight irregularities in the SBM-i results, which is due to the aforementioned mesh sensitivity of the SBM implementations.
<figure><img src="image_11.png" /><figcaption>Figure 9: h-convergence study for the consistent Robin formulations (Aubin and Nitsche) with different $\varepsilon = \{10^{-1}, 10^{0}, 10^{1}\}$ and $P = \{2, 6, 10\}$. The analysis is with CBM and SBM-i.</figcaption></figure>

### 7.3.3 A mixed Dirichlet-Neumann problem using only Robin boundary conditions

To demonstrate the implemented *all-in-one* strategy, we perform a *p*-convergence study on the mesh layout given in Figure 10. Here, a problem with outer Dirichlet and inner Neumann conditions is considered by setting $\varepsilon$ properly in the Robin formulation. In this case, $\varepsilon = 10^{-10}$ on the outer square and $\varepsilon = 10^{10}$ on the inner circle.

In Figure 11, we perform the usual *p*-convergence analysis on three different meshes. Here, spectral convergent behavior for all computations is observed. Moreover, *q* ruined convergence rate for very high orders using SBM-e – due to ill-conditioning – is seen as well. The SBM-i provides similar-sized errors independent of choosing Aubin or Nitsche.
<figure><img src="image_12.png" /><figcaption>Figure 10: Mesh layout with elements of size $l_c = 1/8$ of a (1.5 × 1.5) square domain with a circular hole of radius $R = 0.375$ at the center $(x_0, y_0) = (1, 1)$ is considered. Background mesh domain is of lengths $L_x = L_y = 2$. Robin boundary conditions with $\varepsilon = 10^{-10}$ (Dirichlet mimic) on square and $\varepsilon = 10^{10}$ (Neumann mimic) on circle.</figcaption></figure>

### 7.3.4 Confirmation of asymptotic preserving properties

As a last numerical experiment, we seek to confirm the asymptotic preserving properties of the implementation. Recall the solution perturbations in Proposition 2.1 and Proposition 2.2 for the full Robin solution, $u$. For completeness, in the Dirichlet limit, we have the following expansion (only including up to second-order terms) as

$$u = u_0 + \varepsilon u_1 + \varepsilon^2 u_2 + \mathcal{O}(\varepsilon^3).$$

In this perturbed setting, we solve the full Robin problem for $u$ (for some arbitrary $\varepsilon$) and the zeroth order Dirichlet problem for $u_0$ on a mesh given some approximation order, $P$. Then, $u_1$ and $u_2$ are solved recursively as given in Proposition 2.1. From these solutions, we compute $\|u - u_0\|_{L^1}$, $\|u - u_0 - \varepsilon u_1\|_{L^1}$, and $\|u - u_0 - \varepsilon u_1 - \varepsilon^2 u_2\|_{L^1}$ for different values of $\varepsilon$. The results can be seen in Figure 12, for $P = \{1, ..., 10\}$ on a mesh with $l_c = 1/8$. The convergence studies are supported by contour plots of $u - u_0$ versus $\varepsilon u_1$ and $u - u_0 - \varepsilon u_1$ versus $\varepsilon^2 u_2$ for different

<figure><img src="image_13.png" /><figcaption>Figure 11: $p$-convergence study on mesh layout presented in Figure 10 with $l_c = \{1/4, 1/8, 1/16\}$. For both Aubin and Nitsche with SBM-e and SBM-i.</figcaption></figure>
values of $\epsilon$: $\epsilon = 10^{-2}$ in Figure 13, $\epsilon = 10^{-4}$ in Figure 14, and $\epsilon = 10^{-6}$ in Figure 15. For all contour plots, $P = 10$, where the solution is interpolated onto the true circle for the purpose of visualization.

<figure><img src="image_14.png" /><figcaption>Figure 12: p-convergence on a mesh with $l_c = 1/8$ and variable $\epsilon$ in $\|u - u_0\|_{L^1}$, $\|u - u_0 - \epsilon u_1\|_{L^1}$, and $\|u - u_0 - \epsilon u_1 - \epsilon^2 u_2\|_{L^1}$.</figcaption></figure>

<figure><img src="image_15.png" /><figcaption>Figure 13: Contour plot of $u - u_0$ versus $\epsilon u_1$ and $u - u_0 - \epsilon u_1$ versus $\epsilon^2 u_2$ for $\epsilon = 10^{-2}$ on ($l_c, P) = (1/8, 10$).</figcaption></figure>

From Figure 12, we see how the zeroth order error, $\|u - u_0\|_{L^1}$, converges linearly with $\epsilon$. Moreover, $\|u - u_0 - \epsilon u_1\|_{L^1}$ and $\|u - u_0 - \epsilon u_1 - \epsilon^2 u_2\|_{L^1}$, can be observed to converge as expected second and third order with $\epsilon$, respectively. That is until lower order errors start to dominate, e.g., see for $P = 10$ in the figures. The convergence rate attends spectral behavior until this point. From the contour plots, going through Figure 13 to 15, it is evident how the minimum and maximum values – in general – are attained at the boundary of the domain, whereas it is zero in the center. From
the figures, $u - u_0$ and $\varepsilon u_1$ look to have more and more sufficient agreement as $\varepsilon$ decreases (minor difference for $\varepsilon = 10^{-2}$ and no visible difference for $\varepsilon = 10^{-4}$ and $\varepsilon = 10^{-6}$). Moreover, when decreasing $\varepsilon$ the difference between $u - u_0 - \varepsilon u_1$ and $\varepsilon^2 u_2$ becomes more noticeable. At $\varepsilon = 10^{-6}$, $\varepsilon^2 u_2$ is practically machine precision.

<figure><img src="image_16.png" /><figcaption>Figure 14: Contour plot of $u - u_0$ versus $\varepsilon u_1$ and $u - u_0 - \varepsilon u_1$ versus $\varepsilon^2 u_2$ for $\varepsilon = 10^{-4}$ on $(l_c, P) = (1/8, 10)$.</figcaption></figure>

<figure><img src="image_17.png" /><figcaption>Figure 15: Contour plot of $u - u_0$ versus $\varepsilon u_1$ and $u - u_0 - \varepsilon u_1$ versus $\varepsilon^2 u_2$ for $\varepsilon = 10^{-6}$ on $(l_c, P) = (1/8, 10)$.</figcaption></figure>
# 8 Conclusion

In the paper, we have presented a new high-order spectral element-based numerical model for solving the two-dimensional Poisson problem on an unfitted/immersed/embedded domain. The model utilizes the general idea from the shifted boundary method; however, the cumbersome evaluation of high-order Taylor series expansions is avoided by exploring a polynomial correction. Boundary conditions – Dirichlet, Neumann, and Robin – are enforced weakly in both conformal and unfitted settings. In the latter, we have presented a consistent asymptotic preserving Robin formulation for all values of the Robin parameter, $\varepsilon$. Moreover, we have studied the influence of how to construct the surrogate domain. Ultimately showing that interpolating approaches improve convergence rates and decrease conditioning numbers significantly. Various variational formulations – Aubin and Nitsche type – have been implemented concluding that as long as the basis functions are being interpolated, the differences between the two are minor.

# Acknowledgments

This work contributes to the activities in the PhD research project "New advanced simulation techniques for wave energy converters" hosted at the Department of Applied Mathematics and Computer Science, Technical University of Denmark, Denmark. Work started during a visit for JV at INRIA at the University of Bordeaux.
# A Verification on aligned meshes for Dirichlet problems

## A.1 Aubin's method

<figure><img src="image_18.png" /><figcaption>Figure 16: Verification of the Dirichlet Aubin formulation. From the top row and down: CBM, SBM-e, SBM-ei, and SBM-i. Left column is $h$-convergence in $\|u_h - u\|_{L^1}$ and right column is the associated matrix conditioning, $\kappa(\mathcal{A})$. Analysis is for $P = \{1, ..., 5\}$ on aligned meshes.</figcaption></figure>
A.2 Nitsche's method

<figure><img src="image_19.png" /><figcaption>Figure 17: Verification of the Dirichlet Nitsche formulation. From the top row and down: CBM, SBM-e, SBM-ei, and SBM-i. Left column is $h$-convergence in $\|u_h - u\|_{L^1}$ and right column is the associated matrix conditioning, $\kappa(\mathcal{A})$. Analysis is for $P = \{1, \dots, 5\}$ on well-behaved meshes.</figcaption></figure>
## B Assessment on unstructured meshes for Robin boundary conditions

<figure><img src="image_20.png" /><figcaption>Figure 18: Median (top row) and variance (bottom row) of 100 discrete solutions with the setup as in Figure 5 for the consistent Robin problem with $\epsilon = 1$. Quantities are for $||Au - b||_{L^1}$ and $\kappa(A)$ for Aubin and Nitsche formulation with SBM-e, SBM-ei, and SBM-i.</figcaption></figure>

## References

[1] J. Li, S. Rosencrans, X. Wang, and K. Zhang. Asymptotic analysis of a Dirichlet problem for the heat equation on a coated body. *Proc.Amer.Math.Soc.*, 137(5):1711–1721, 2008.

[2] H. Li and X. Wang. Effective boundary conditions for the heat equation with interior inclusion. *Commun.Math.Res.*, 36(3):272–295, 2020.

[3] Yves Achdou, O. Pironneau, and F. Valentin. Effective boundary conditions for laminar flows over periodic rough boundaries. *J. Comput. Phys.*, 147(1):187–218, 1998.

[4] X. Antoine and H. Barucq. Approximation by generalized impedance boundary conditions of a transmissions problem in acoustic scattering. *ESAIM: M2AN*, 39:1041–1059, 2005.

[5] J. Dias and V. Peron. Equivalent Robin boundary conditions for acoustic and elastic media. *Math Models Methods Appl Sci.*, 28(8):1531–1566, 2016.

[6] Finnur Pind, Cheol-Ho Jeong, Jan S. Hesthaven, Allan P. Engsig-Karup, and Jakob Strømann-Andersen. A phenomenological extended-reaction boundary model for time-domain wave-based acoustic simulations under sparse reflection conditions using a wave splitting method. *Applied Acoustics*, 172:107596, 2021.

[7] H. Ammari and S. He. Effective imedence boundary conditions for an inhomogenous thin layer on a curved metallig surface. *IEEE Trans. Antennas Propag.*, 46(5):710–715, 1998.
[8] N. Bartoli and A. Bendali. Robust and high order effusive boundary conditions for perfectly conducting scatterers coated by a thin dielectonic layer. *IMA J. Appl. Math.*, 67:479–508, 2002.

[9] H. Haddar, P. Joly, and H.-M. Nguyen. Generalized impedance boundary conditions for scattering by strongly absorbing obstacles: the scalar case. *Math. Mod. Meth. Appl. Sci.*, 15(8):1273–1300, 2005.

[10] B. Aslanurek, H. Haddar, and H. Sahinturk. Generalized impedance boundary conditions for thin dielectric coatings with variable thickness. *Wave Motion*, 48:681–700, 2011.

[11] Ronan Perrussel and Clair Poignard. Asymptotic expansion of steady-state potential in a high contrast medium with a thin resistive layer. *Appl. Math. Comput.*, 221:48–65, 2013.

[12] A. Auvray and G. Vial. Asymptotic expansions and effective boundary conditions: a short review for smooth and non-smooth geometries with thin layers. *ESAIM: Proceedings and Surveys*, 61:38–54, 2018.

[13] Luca Gerardo-Giorda, Fabio Nobile, and Christian Vergara. Analysis and optimization of robin-robin partitioned procedures in fluid-structure interaction problems. *SIAM J. Numer. Anal.*, 48(6):2091–2116, 2010.

[14] Santiago Badia, Fabio Nobile, and Christian Vergara. Fluid-structure partitioned procedures based on robin transmission conditions. *J. Comput. Phys.*, 227(14):7027–7051, 2008.

[15] Tagliabue, Anna, Dedè, Luca, and Quarteroni, Alfio. Nitsche's method for parabolic partial differential equations with mixed time varying boundary conditions. *ESAIM: M2AN*, 50(2):541–563, 2016.

[16] Heinz-Otto Kreiss and Joseph Oliger. Comparison of accurate methods for the integration of hyperbolic equations. *Tellus*, 24(3):199–215, 1972.

[17] Jan S. Hesthaven, Sigal Gottlieb, and David Gottlieb. *Spectral Methods for Time-Dependent Problems*. Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press, 2007.

[18] Hui Xu, Chris Cantwell, Carlos Monteserin, Claes Eskilsson, Allan Peter Engsig-Karup, and Spencer Sherwin. Spectral/hp element methods: Recent developments, applications, and perspectives. *J. Hydrodynam. B*, 30, 02 2018.

[19] Ivo Babuška and Manil Suri. The p- and h-p versions of the finite element method, an overview. *Comput. Methods. Appl. Mech. Eng.*, 80(1):5–26, 1990.

[20] Anthony T Patera. A spectral element method for fluid dynamics: Laminar flow in a channel expansion. *J. Comput. Phys.*, 54(3):468–488, 1984.

[21] E. M. Rønquist and A. T. Patera. Spectral element multigrid. i. formulation and numerical results. *J. Sci. Comput.*, 2:389–406, 1987.

[22] A. P. Engsig-Karup and W.L. Laskowski. An efficient *p*-multigrid spectral element model for fully nonlinear water waves and fixed bodies. *Int J Numer Methods Fluids*, 93(9):2823–2841, 2021.

[23] Charles S Peskin. Flow patterns around heart valves: A numerical method. *J. Comput. Phys.*, 10(2):252–271, 1972.

[24] Roberto Verzicco. Immersed boundary methods: Historical perspective and future outlook. *Annu. Rev. Fluid Mech.*, 55(1):129–155, 2023.

[25] Ramakrishnan Thirumalaisamy, Neelesh A. Patankar, and Amneet Pal Singh Bhalla. Handling neumann and robin boundary conditions in a fictitious domain volume penalization framework. *J. Comput. Phys.*, 448:110726, 2022.

[26] Isabelle Ramière, Philippe Angot, and Michel Belliard. A fictitious domain approach with spread interface for elliptic problems with general boundary conditions. *Comput. Methods. Appl. Mech. Eng.*, 196(4):766–781, 2007.

[27] Bouchra Bensiali, Guillaume Chiavassa, and Jacques Liandrat. Penalization of robin boundary conditions. *Appl. Numer Math*, 96:134–152, 2015.

[28] Rajat Mittal and Gianluca Iaccarino. Immersed boundary methods. *Annu. Rev. Fluid Mech.*, 37:239–261, 1 2005.

[29] T Strouboulis, I Babuška, and K Copps. The design and analysis of the generalized finite element method. *Comput. Methods. Appl. Mech. Eng.*, 181:43–69, 2000.

[30] Nicolas Moës, John Dolbow, and Ted Belytschko. A finite element method for crack growth without remeshing. *Int. J. Numer. Methods. Eng.*, 46:131–150, 9 1999.

[31] Anita Hansbo and Peter Hansbo. An unfitted finite element method, based on Nitsche's method, for elliptic interface problems. *Comput. Methods. Appl. Mech. Eng.*, 191(47):5537–5552, 2002.

[32] Erik Burman, Susanne Claus, Peter Hansbo, Mats G. Larson, and André Massing. Cutfem: Discretizing geometry and partial differential equations. *Int. J. Numer. Methods. Eng.*, 104(7):472–501, 2015.
[33] Santiago Badia, Francesc Verdugo, and Alberto F. Martín. The aggregated unfitted finite element method for elliptic problems. *Comput. Methods. Appl. Mech. Eng.*, 336:533–553, 2018.

[34] Christoph Lehrenfeld. High order unfitted finite element methods on level set domains using isoparametric mappings. *Comput. Methods. Appl. Mech. Eng.*, 300:716–733, 2016.

[35] R. Oyarzúa, M. Solano, and P. Zúniga. A high order mixed-fem for diffusion problems on curved domains. *J.Sci.Comput.*, 79:49–78, 2019.

[36] Yuanming Xiao, Jinchao Xu, and Fei Wang. High-order extended finite element methods for solving interface problems. *Comput. Methods. Appl. Mech. Eng.*, 364:112964, 2020.

[37] Alexei Lozinski. Cutfem without cutting the mesh cells: A new way to impose dirichlet and neumann boundary conditions on unfitted meshes. *Comput. Methods. Appl. Mech. Eng.*, 356:75–100, 2019.

[38] Michel Duprez, Vanessa Lleras, and Alexei Lozinski. A new $\phi$-fem approach for problems with natural boundary conditions. *Numer. Methods. Partial. Differ. Equ.*, 39(1):281–303, 2023.

[39] Alex Main and Guglielmo Scovazzi. The shifted boundary method for embedded domain computations. part i: Poisson and stokes problems. *J. Comput. Phys.*, 372:972–995, 2018.

[40] Alex Main and Guglielmo Scovazzi. The shifted boundary method for embedded domain computations. part ii: Linear advection-diffusion and incompressible navier-stokes equations. *J. Comput. Phys.*, 372:996–1026, 2018.

[41] 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 flows. *J. Comput. Phys.*, 369:45–79, 2018.

[42] Mirco Ciallella, Elena Gaburro, Marco Lorini, and Mario Ricchiuto. Shifted boundary polynomial corrections for compressible flows: high order on curved domains using linear meshes. *Appl. Math. Comput.*, 441:127698, 2023.

[43] Nabil M. Atallah, Ketan Mittal, Guglielmo Scovazzi, and Vladimir Z. Tomov. A high-order shifted interface method for lagrangian shock hydrodynamics, 2023.

[44] Kangan Li, John G. Michopoulos, Athanasios Iliopoulos, John C. Steuben, and Guglielmo Scovazzi. Complex-geometry simulations of transient thermoelasticity with the shifted boundary method. *Comput. Methods. Appl. Mech. Eng.*, 418:116461, 2024.

[45] Kangan Li, Nabil M Atallah, G Alex Main, and Guglielmo Scovazzi. The shifted interface method: A flexible approach to embedded interface computations. *Int. J. Numer. Methods. Eng.*, 121:492–518, 2020.

[46] Oriol Colomés, Alex Main, Léo Nouveau, and Guglielmo Scovazzi. A weighted shifted boundary method for free surface flow problems. *J. Comput. Phys.*, 424:109837, 2021.

[47] T. Carlier, L. Nouveau, H. Beaugendre, M. Colin, and M. Ricchiuto. An enriched shifted boundary method to account for moving fronts. *J. Comput. Phys.*, 489:112295, 2023.

[48] Nabil M Atallah, Claudio Canuto, and Guglielmo Scovazzi. The shifted boundary method for solid mechanics. *Int. J. Numer. Methods. Eng.*, 122:5935–5970, 2021.

[49] Kangan Li, Nabil M Atallah, Antonio Rodríguez-Ferran, Dakshina M Valiveti, and Guglielmo Scovazzi. The shifted fracture method. *Int. J. Numer. Methods. Eng.*, 122:6641–6679, 2021.

[50] Leo Nouveau, Mario Ricchiuto, and Guglielmo Scovazzi. High-order gradients with the shifted boundary method: an embedded enriched mixed formulation for elliptic pdes. *J. Comput. Phys.*, 398:108898, 2019.

[51] Nabil M Atallah, Claudio Canuto, and Guglielmo Scovazzi. The high-order shifted boundary method and its analysis. *Comput. Methods. Appl. Mech. Eng.*, 394:114885, 2022.

[52] J. Haydel Collins, Alexei Lozinski, and Guglielmo Scovazzi. A penalty-free shifted boundary method of arbitrary order. *Comput. Methods. Appl. Mech. Eng.*, page 116301, 2023.

[53] Cheng-Hau Yang, Kumar Saurabh, Guglielmo Scovazzi, Claudio Canuto, Adarsh Krishnamurthy, and Baskar Ganapathysubramanian. Optimal surrogate boundary selection and scalability studies for the shifted boundary method on octree meshes, 2023.

[54] Santiago Badia, Eric Neiva, and Francesc Verdugo. Robust high-order unfitted finite elements by interpolation-based discrete extension. *Comput. Math. with Appl.*, 127:105–126, 2022.

[55] Joachim Nitsche. Über ein variationsprinzip zur lösung von dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind. *Abhandlungen aus dem mathematischen Seminar der Universität Hamburg*, 36:9–15, 1971.
[56] Joseph Benzaken, John A Evans, and Rasmus Tamstorf. Constructing nitsche's method for variational problems. *arXiv preprint arXiv:2203.02603*, 2022.

[57] Mika Juntunen and Rolf Stenberg. Nitsche's method for general boundary conditions. *Math. Comput.*, 78(267): 1353–1374, 2009.

[58] M. Winter, B. Schott, A. Massing, and W.A. Wall. A nitsche cut finite element method for the oseen problem with general navier boundary conditions. *Comput. Methods. Appl. Mech. Eng.*, 330:220–252, 2018.

[59] Jean-Pierre Aubin. Approximation des problèmes aux limites non homogènes pour des opérateurs non linéaires. *J. Math. Anal. Appl.*, 30(3):510–521, 1970.

[60] Douglas N. Arnold, Franco Brezzi, Bernardo Cockburn, and L. Donatella Marini. Unified analysis of discontinuous galerkin methods for elliptic problems. *SIAM J. Numer. Anal.*, 39(5):1749–1779, 2002.

[61] J.D. Murray. *Asymptotic analysis*. Springer New York, NY, 1984.

[62] Radyadour Kh. Zeytounian. *Asymptotic Modelling of Fluid Flow Phenomena*. Springer Dordrecht, 2002.

[63] Stanley Osher and Ronald P. Fedkiw. Level set methods: An overview and some recent results. *J. Comput. Phys.*, 169(2):463–502, 2001.

[64] Roland Glowinski, Tsorng-Whay Pan, and Jacques Periaux. A fictitious domain method for dirichlet problem and applications. *Comput. Methods. Appl. Mech. Eng.*, 111(3):283–303, 1994.

[65] Erik Burman, Peter Hansbo, Mats G. Larson, André Massing, and Sara Zahedi. Full gradient stabilized cut finite element methods for surface partial differential equations. *Comput. Methods. Appl. Mech. Eng.*, 310:278–296, 2016.

[66] Mats G Larson and Sara Zahedi. Stabilization of high order cut finite element methods on surfaces. *SIAM J. Numer. Anal.*, 40(3):1702–1745, 04 2019.

[67] Christoph Lehrenfeld, Maxim A. Olshanskii, and Xianmin Xu. A stabilized trace finite element method for partial differential equations on evolving surfaces. *SIAM J. Numer. Anal.*, 56(3):1643–1672, 2018.

[68] Benoit Fabrèges, Loïc Gouarin, and Bertrand Maury. A smooth extension method. *Comptes Rendus. Mathématique*, 351(9-10):361–366, 2013.

[69] George Karniadakis and Spencer Sherwin. *Spectral/hp element methods for computational fluid dynamics*. OUP Oxford, 2005.

[70] Jan S Hesthaven and Tim Warburton. *Nodal discontinuous Galerkin methods: algorithms, analysis, and applications*. Springer Science & Business Media, 2007.

[71] Tim Warburton. An explicit construction of interpolation nodes on the simplex. *J. Eng. Math.*, 56:247–262, 2006.

[72] Frits de Prenter, Christoph Lehrenfeld, and André Massing. A note on the stability parameter in nitsche's method for unfitted boundary value problems. *Comput. Math. with Appl.*, 75(12):4322–4336, 2018.

[73] A. Ern and J.-L. Guermond. *Theory and practice of finite elements*, volume 159 of *Appl. Math. Sci*. Springer, 2004.

[74] A. Ern and J.-L. Guermond. *Finite Elements II: Galerkin Approximation, Elliptic and MixedPDEs*. Texts in Applied Mathematics. Springer, 2021.

[75] Nabil M Atallah, Claudio Canuto, and Guglielmo Scovazzi. The second-generation shifted boundary method and its numerical analysis. *Comput. Methods. Appl. Mech. Eng.*, 372:113341, 2020.

[76] Erik Burman. A penalty-free nonsymmetric nitsche-type method for the weak imposition of boundary conditions. *SIAM J. Numer. Anal.*, 50(4):1959–1981, 2012.

[77] Richard Pasquetti and Francesca Rapetti. Spectral element methods on unstructured meshes: which interpolation points? *Numer. Algorithms*, 55:349–366, 2010.

[78] Qi Chen and Ivo Babuška. Approximate optimal points for polynomial interpolation of real functions in an interval and in a triangle. *Comput. Methods. Appl. Mech. Eng.*, 128(3):405–417, 1995.

[79] Christophe Geuzaine and Jean-François Remacle. Gmsh: A 3-d finite element mesh generator with built-in pre- and post-processing facilities. *Int. J. Numer. Methods. Eng.*, 79(11):1309–1331, 2009.

[80] Michel Duprez, Vanessa Lleras, and Alexei Lozinski. A new $\phi$-fem approach for problems with natural boundary conditions. *Numer. Methods. Partial. Differ. Equ.*, 39(1):281–303, 2023.