# A nodal ghost method based on variational formulation and regular square
grid for elliptic problems on arbitrary domains in two space dimensions

Clarissa Astuto  King Abdullah University of Science and Technology (KAUST),
4700, Thuwal, Saudi Arabia  Daniele Boffi  King Abdullah University of Science
and Technology (KAUST), 4700, Thuwal, Saudi Arabia  Dipartimento di Matematica
”F. Casorati”, University of Pavia, 27100, Pavia, Italy  Giovanni Russo
Department of Mathematics and Computer Science, University of Catania, 95125,
Catania, Italy  Umberto Zerbinati  University of Oxford, OX2 6GG, Oxford,
England

###### Abstract

This paper focuses on the numerical solution of elliptic partial differential
equations (PDEs) with Dirichlet and mixed boundary conditions, specifically
addressing the challenges arising from irregular domains. Both finite element
method (FEM) and finite difference method (FDM), face difficulties in dealing
with arbitrary domains. The paper introduces a novel nodal symmetric ghost
finite element method approach, which combines the advantages of FEM and FDM.
The method employs bilinear finite elements on a structured mesh, and provides
a detailed implementation description. A rigorous a priori convergence rate
analysis is also presented. The convergence rates are validated with many
numerical experiments, in both one and two space dimensions.

##  1 Introduction

Elliptic partial differential equations (PDEs) are used to describe a large
variety of physical phenomena with applications ranging from modeling the
diffusion of a certain chemical species in a solvent to determining the
electric field corresponding to a given charge density. In the interest of
keeping our presentation as streamlined as possible we will only focus our
attention on the simplest elliptic PDE, i.e. the Poisson equation with mixed
boundary conditions,

$-\Delta u=f\;{\rm in}\;\Omega,\;u=g_{D}\;{\rm
on}\;\Gamma_{D}\subset\partial\Omega\;{\rm and}\;\frac{\partial u}{\partial
n}=g_{N}\;{\rm on}\;\Gamma_{N}\subset\partial\Omega.$ (1)

The previous equation is defined on a domain $\Omega\subset
R\subseteq\mathbb{R}^{d}$, where $R$ is a rectangular region, and $d=1,2$ is
the number of space dimensions.

When solving (1) in two space dimensions, a significant difficulty comes from
the arbitrary shape of $\Omega$. Indeed, while for simple geometries, such as
a square or a circle, it is possible to find an analytical solution in terms
of trigonometric and Bessel functions, for more complex geometries this is not
possible and one has to resort to numerical methods.

Among the several schemes developed to solve (1) on arbitrary domains, two
methods of great importance are the finite element method (FEM) and the finite
difference method (FDM).

The body fitted FEM is a very general method that relies on a tessellation of
the domain $\Omega$, the construction of “good” tessellation is far from being
a trivial task, [25, 26, 32, 42]. Since the generation of a mesh for $\Omega$
is such a delicate task, the use of descriptions that are closely related to
the geometric data acquisition can dramatically reduce the cost of
preprocessing an acquired geometry into representations suitable for the FEM.

For example, the generalization of standard finite elements, known as
isogeometrical analysis, uses a choice of basis function compatible with the
description of the geometry adopted by most CAD software, [33, 24]. Another
extremely powerful approaches, known as the CutFEM, are finite element
discretizations where the boundary of $\Omega$ is represented on a background
mesh, which can be of low-quality, using a level set function. In the CutFEM
the background mesh is also used to represent the numerical solution of (1),
[14, 34, 35, 16]. Since the CutFEM does not require the use of a mesh fitted
to the body over which we are solving the Poisson equation, it belongs to the
realm of the “unfitted” finite elements methods.

An hybrid approach between fitted and unfitted finite elements methods is the
Fictitious Domain method, first presented in [29]. In this method a body
fitted mesh is only constructed for boundary $\partial\Omega$ and a background
mesh of the computational domain is used to represent the numerical solution.
In the Fictitious Domain method the boundary condition are imposed via
distributed Lagrange multipliers, similarly to [9] and [10]. A comparative
study between CutFEM and the Lagrange multiplier method previously discussed
can be found in [8].

The CutFEM is just one of the methods that makes use of level set to describe
the geometry of the domain. Level set in general will allow to describe the
geometry of the domain in a very flexible way, [39, 45, 40].

Not only the FEM can be used in combination with level set techniques to
describe the geometry of the domain, but there is also a vast literature of
unfitted finite difference schemes [28, 17, 44, 19, 20], just to mention a
few. A common feature of most unfitted FDMs previously mentioned is the
addition of ghost points, i.e. points that are not part of the domain $\Omega$
but are used to impose the boundary conditions and retrieve optimal order of
convergence. We can find many applications in [21, 4, 1, 2, 3].

While the unfitted FDM, compared to unfitted FEM, might be easier to implement
and understand for a novice user, these finite difference schemes have the
disadvantage of preserving less structural properties of the continuous
problem, such as the symmetry and the positive definiteness of the stiffness
matrix. In the present paper we aim to present ideas from the unfitted FEM and
FDM literature and introduce the reader to a nodal ghost finite element method
for elliptic problems on irregular domains with mixed boundary conditions.

The choice of bilinear finite elements on a structured mesh yield a similar
stencil to the one typical of FDM, allowing a much more natural and straight-
forward implementation for a novice user.

To achieve this goal, the paper is structured as follows: we begin deriving
the variational formulation of (1) in the one dimensional case and describe
our nodal symmetric ghost point method for this simple case. We then proceed
to extend our nodal symmetric ghost point method to the two dimensional case.
In Section 5, we provide a detailed description of the implementation aspects
of the scheme we propose. In particular, we discuss the exact quadrature of
the boundary integral, which is one of the main novelty of our work. In
Section 6, we prove a priori convergence rate for the proposed scheme.
Furthermore, to the best our knowledge we give for the first time a rigorous
treatment of the snapping back to grid perturbation. Lastly, we validate the
method and the a priori estimated we proved via extensive numerical
experiments in both one and two space dimensions.

We hope that our treatment of the proposed nodal symmetric ghost point method,
which starts from the most basic concepts and build up to more complicated
scenario, will facilitate the novice reader understating of our scheme.

##  2 The continuous one dimensional problem

Let us start our discussion from the Poisson equation with homogeneous
Dirichlet conditions at the boundary defined in a domain $\Omega=[a,b]\subset
R=[0,1]$:

$-u^{\prime\prime}(x)=f(x),\quad x\in[a,b],\quad u(a)=0,\>u(b)=0.$ (2)

Multiplying the above equation by a test function $v\in H^{1}_{0}([a,b])$,
provides the following variational formulation of the problem:

$\int_{a}^{b}\left[u^{\prime}(x)v^{\prime}(x)-f(x)v(x)\right]\,dx=u^{\prime}(b)v(b)-u^{\prime}(a)v(a),\quad\forall
v\in H^{1}_{0}([a,b])$ (3)

where $H^{1}_{0}([a,b])$ is the standard Sobolev space with homogeneous trace
at the boundary [12]. Thanks to this condition, the right hand side of the
equation vanishes, leading to the formulation:

$\int_{a}^{b}\left[u^{\prime}(x)v^{\prime}(x)-f(x)v(x)\right]\,dx=0,\quad\forall
v\in H^{1}_{0}([a,b])$ (4)

The variational formulation can be obtained setting to zero the first
variation of the functional

$J_{D0}[u]=\int_{a}^{b}\left[\frac{1}{2}(u^{\prime}(x))^{2}\,dx-f(x)u(x)\right]\,dx,$
(5)

as can be immediately verified by performing the first variation of the
functional by replacing $u$ by $u+\varepsilon v$, with $u$ and $v$ in
$H^{1}_{0}([a,b])$, and imposing

$\left.\frac{d}{d\varepsilon}J_{D0}[u+\varepsilon
v]\right|_{\varepsilon=0}=0\quad\forall v\in H^{1}_{0}([a,b]).$ (6)

In particular we notice that since, $J_{D0}$ is a convex functional, looking
for functions for which the first variation of $J_{D0}$ vanishes is equivalent
to minimizing the functional.

A different way to approximately impose zero Dirichlet conditions on the
boundary is to penalise any discrepancy from zero. This is obtained by adding
a penalisation term to the functional:

$J_{D0\lambda}[u]=\int_{a}^{b}\left[\frac{1}{2}(u^{\prime}(x))^{2}\,dx-f(x)u(x)\right]\,dx+\frac{\lambda}{2}[u(a)^{2}+u(b)^{2}].$
(7)

Notice that the energy functional $J_{D0\lambda}$ corresponds to (2) where the
Dirichlet boundary conditions are replaced by Robin boundary conditions. It is
a well-established fact that as the penalisation grows to infinity it will
force zero Dirichlet condition as $\lambda\to\infty$ [22, 23]. When
integrating by parts, we obtain the following variational formulation of the
problem:

$\displaystyle\int_{a}^{b}\left[u^{\prime}(x)v^{\prime}(x)-f(x)v(x)\right]\,dx$
$\displaystyle+[\lambda u(a)+u^{\prime}(a)]v(a)$ (8) $\displaystyle+[\lambda
u(b)-u^{\prime}(b)]v(b)=0,\quad\forall v\in H^{1}([a,b])$

where now we do not need to impose that $v$ is zero on the boundary. Notice
that the variational formulation (8) is the starting point for the penalty
finite element proposed in [5, 6] and analyzed in [36, 7].

If we want to impose non homogeneous Dirichlet conditions, say we impose
$u(a)=u_{a}$, $u(b)=u_{b}$, we have to consider the functional

$J_{D\lambda}[u]=\int_{a}^{b}\left[\frac{1}{2}(u^{\prime}(x))^{2}-f(x)u(x)\right]\,dx+\frac{\lambda}{2}[(u(a)-u_{a})^{2}+(u(b)-u_{b})^{2}]$
(9)

which leads to the following variational formulation:

$\begin{split}\int_{a}^{b}\left[u^{\prime}(x)v^{\prime}(x)\right.&\left.-f(x)v(x)\right]\,dx+[\lambda(u(a)-u_{a})+u^{\prime}(a)]v(a)\\\
&+[\lambda(u(b)-u_{b})-u^{\prime}(b)]v(b)=0,\quad\forall v\in
H^{1}([a,b])\end{split}$ (10)

The variational formulation (10) is non symmetric, because of the term
$\displaystyle\left.(u^{\prime}v)\right|_{a}^{b}$. We will now address this
issue introducing an alternative formulation, which is equivalent at a
continuous level and is obtained as follows. We add and subtract the quantity
$\displaystyle\left.(uv^{\prime})\right|_{a}^{b}$:

$\int_{a}^{b}u^{\prime}(x)v^{\prime}(x)\,dx-\left.(u^{\prime}v)\right|_{a}^{b}-\left.(uv^{\prime})\right|_{a}^{b}=\int_{a}^{b}f(x)v(x)\,dx-
u_{b}v^{\prime}(b)+u_{a}v^{\prime}(a)$ (11)

where we made use of the Dirichlet boundary conditions in $a$ and $b$ on the
right-hand side.

The coercivity of the operator may be improved by adding a penalisation term,
therefore for continuous formulation can be written as

$\displaystyle\int_{a}^{b}u^{\prime}(x)v^{\prime}(x)\,dx-\left.(u^{\prime}v)\right|_{a}^{b}-\left.(uv^{\prime})\right|_{a}^{b}+$
$\displaystyle\lambda(u(a)-u_{a})v(a)+\lambda(u(b)-u_{b})v(b)$
$\displaystyle=\int_{a}^{b}f(x)v(x)\,dx-u_{b}v^{\prime}(b)+u_{a}v^{\prime}(a)$
(12)

and bringing to the right-hand side the known term, we have

$\displaystyle\int_{a}^{b}u^{\prime}(x)v^{\prime}(x)\,dx-\left.(u^{\prime}v)\right|_{a}^{b}-\left.(uv^{\prime})\right|_{a}^{b}+\lambda
u(a)v(a)+\lambda u(b)v(b)$ $\displaystyle=$
$\displaystyle\int_{a}^{b}f(x)v(x)\,dx-
u_{b}v^{\prime}(b)+u_{a}v^{\prime}(a)+\lambda u_{a}v(a)+\lambda u_{b}v(b)$
(13)

Is trivial to obtain the analogous variational formulation for the Neumann
problem, in fact if we consider

$-u^{\prime\prime}(x)=f(x),\quad x\in[a,b]\qquad u^{\prime}(a)=g_{a},\quad
u^{\prime}(b)=g_{b},$ (14)

we can multiply by a test function $v\in H^{1}[a,b]$, integrate over $[a,b]$
and integrate by parts to see that the solution of (14)$is$ $u\in H^{1}[a,b]$
such that

$\int_{a}^{b}u^{\prime}(x)v^{\prime}(x)\,dx=\int_{a}^{b}f(x)v(x)\,dx+g_{b}v(b)-g_{b}av(a)\quad\forall
v\in H^{1}([a,b]).$ (15)

The only aspects we need to be careful when considering a variational
formulation for the Neumann problem is that a compatibility condition has to
be verified by $f\in L^{2}([a,b])$ and we need to choose a unique solution in
the kernel of the gradient, imposing a null-average condition, i.e.

$\int_{a}^{b}f\,dx=g_{b}-g_{a},\qquad\int_{a}^{b}u\,dx=0.$ (16)

Let us now consider the following mixed problem:

$-u^{\prime\prime}(x)=f(x),\quad x\in[a,b]\qquad u(a)=0,\quad
u^{\prime}(b)=g_{b}$ (17)

We shall look for solutions of (17) in $H^{1}([a,b])$, respecting the boundary
conditions. Multiplying by a test function

$H^{1}_{a}([a,b]):=\left\\{v\in H^{1}([a,b])\;:\;v(a)=0\right\\},$ (18)

and integrating by parts in $[a,b]$ one obtains:

$\int_{a}^{b}u^{\prime}(x)v^{\prime}(x)\,dx-\left.(u^{\prime}v)\right|_{a}^{b}=\int_{a}^{b}f(x)v(x)\,dx.$
(19)

Now, making use of the boundary condition in $b$ and the fact that for any
$v\in H^{1}_{a}([a,b])$ $v(a)=0$, we obtain the following variational
formulation of problem (17):

$\int_{a}^{b}u^{\prime}(x)v^{\prime}(x)\,dx=\int_{a}^{b}f(x)v(x)\,dx+g_{b}v(b)\quad\forall
v\in H^{1}_{a}([a,b]).$ (20)

This equation can be obtained by looking for the function $u\in
H^{1}_{a}([a,b])$ that minimizes the following functional:

$J[u]=\frac{1}{2}\int_{a}^{b}[(u^{\prime}(x))^{2}\,dx-f(x)u(x)]\,dx+g_{b}u(b)$
(21)

Indeed, by equating to zero the first variation of the functional, i.e.
imposing

$\left.\frac{d}{d\varepsilon}J[u+\varepsilon v]\right|_{\varepsilon=0}=0$ (22)

$\forall v\in H^{1}_{a}([a,b])$, gives exactly (20). As in the case of the
Dirichlet problem, we may impose an approximate Dirichlet condition on the
left hand boundary by penalising the functional as follows:

$J[u]=\int_{a}^{b}\left[\frac{1}{2}(u^{\prime}(x))^{2}\,dx-f(x)u(x)\right]\,dx+g_{b}u(b)+\frac{\lambda}{2}(u(a)-u_{a})^{2},\quad
u\in H^{1}([a,b])$

which leads to the following variational formulation of the problem in
$H^{1}([a,b])$:

$\int_{a}^{b}u^{\prime}(x)v^{\prime}(x)\,dx+[\lambda(u(a)-u_{a})+u^{\prime}(a)]v(a)=\int_{a}^{b}f(x)v(x)\,dx+g_{b}v(b).$
(23)

As before, we look for a term that symmetries (23). For this reason, we add
and subtract the term $u(a)v^{\prime}(a)$, making use of the Dirichlet
boundary condition on $a$

$\displaystyle\int_{a}^{b}u^{\prime}(x)v^{\prime}(x)\,dx+[\lambda(u(a)-u_{a})+u^{\prime}(a)]v(a)+u(a)v^{\prime}(a)$
(24) $\displaystyle=\int_{a}^{b}f(x)v(x)\,dx+u_{a}v^{\prime}(a)+g_{b}v(b)$

and bringing to the right-hand side the known terms, we have

$\displaystyle\int_{a}^{b}u^{\prime}(x)v^{\prime}(x)\,dx+\lambda
u(a)v(a)u_{a}+u^{\prime}(a)v(a)+u(a)v^{\prime}(a)$ (25)
$\displaystyle=\int_{a}^{b}f(x)v(x)\,dx+u_{a}v^{\prime}(a)+g_{b}v(b)+\lambda
u_{a}v(a).$

The technique we used to symmetrize the variational problem is known as the
Nitsche technique [38].

##  3 The discrete one-dimensional problem

The region $R=[0,1]$ is discretized by a uniform mesh of grid size $h=1/N$.

Without loss of generality we assume that $a\in[0,h]$, and $b\in[1-h,1]$. If
$a=0$ and $b=1$, standard finite element discretization in the interval
$[0,1]$ can be used.

We denote by $\theta_{1}=1-a/h$ and $\theta_{2}=(1-b)/h$, so that
$\theta_{1},\theta_{2}\in[0,1]$. We denote by $I_{j}=[j-1,j]h$, with
$j=1,\ldots,N$, the $j$-th interval of the subdivision. We consider the space
$V_{h}$ of dimension $N+1$ as the space of the splines $S_{1}$ on the mesh
$\\{I_{j}\\}_{i=1}^{N}$. As basis for the space $V_{h}$ we choose the
functions $\varphi_{i}(x)$, $i=0,\ldots,N,$ defined as

$\varphi_{i}(x)=\max\\{1-|x-x_{i}|/h,0\\},\quad x\in[0,1],\quad i=0,\ldots,N.$
(26)

A schematic representation of the discretization setup is illustrated in
Figure 1.

Figure 1: Discretization setup in 1D: the region $R=[0,1]\supset[a,b]=\Omega$,
and the basis function $\varphi_{j},j=0,\ldots,N$, defined in (26).

A generic function $u_{h}\in V_{h}$ can be expressed using the basis
$\\{\varphi_{i}\\}_{i=1}^{N}$ as follows:

$u_{h}(x)=\sum_{j=1}^{N}u_{j}\varphi_{j}(x)$ (27)

Now we look for a function $u_{h}\in V_{h}$ that solves the variational
formulation of the equation (10) in the case of Dirichlet problem, and (23) in
the case of mixed problem.

We now make use of the conformity of the discrete space, i.e. $V_{h}\subset
H^{1}_{0}([a,b])$, to obtain a discrete variational problem, by considering
(10) with both test and trial functions in $V_{h}$:

$\begin{split}&\int_{a}^{b}u_{h}^{\prime}(x)\varphi_{i}^{\prime}(x)\,dx-\int_{a}^{b}f(x)\varphi_{i}(x)\,dx+\lambda(u_{h}(a)-u_{a})\varphi_{i}(a)+u_{h}^{\prime}(a)\varphi_{i}(a)\\\
&+u_{h}(a)\varphi_{i}^{\prime}(a)+\lambda(u_{h}(b)-u_{b})\varphi_{i}(b)-u_{h}^{\prime}(b)\varphi_{i}(b)-u_{h}(b)\varphi_{i}^{\prime}(b)=0,\quad
i=0,\ldots,N\end{split}$

By replacing $u_{h}$ with its representation (27), we obtain the following
system of equations for the unknowns $u_{i}$, $i=1,\ldots,N$:

$\displaystyle\sum_{j=0}^{N}u_{j}\underbrace{(\varphi_{j}^{\prime},\varphi_{i}^{\prime})_{L^{2}([a,b])}}_{\rm
S}+\sum_{j=0}^{N}u_{j}\underbrace{\left[\varphi_{j}^{\prime}(a)\varphi_{i}(a)+\varphi_{j}(a)\varphi_{i}^{\prime}(a)-\varphi_{j}^{\prime}(b)\varphi_{i}(b)-\varphi_{j}(b)\varphi_{i}^{\prime}(b)\right]}_{\rm
S_{T}}$
$\displaystyle+\lambda\sum_{j=0}^{N}u_{j}\left(\underbrace{\varphi_{j}(a)\varphi_{i}(a)}_{\rm
P_{a}}+\underbrace{\varphi_{j}(b)\varphi_{i}(b)}_{\rm P_{b}}\right)$ (28)
$\displaystyle=\sum_{j=0}^{N}f_{j}\underbrace{(\varphi_{j},\varphi_{i})_{L^{2}([a,b])}}_{\rm
M}\,dx+\sum_{j=0}^{N}{u_{a}}_{j}\underbrace{\varphi_{j}(a)\varphi_{i}^{\prime}(a)}_{\rm
D_{a}}+\sum_{j=0}^{N}{u_{b}}_{j}\underbrace{\varphi_{j}(b)\varphi_{i}^{\prime}(b)}_{\rm
D_{b}}+$ (29)
$\displaystyle+\lambda\left(\sum_{j=0}^{N}{u_{a}}_{j}\underbrace{\varphi_{j}(a)\varphi_{i}(a)}_{\rm
P_{a}}+\sum_{j=0}^{N}{u_{b}}_{j}\underbrace{\varphi_{j}(b)\varphi_{i}(b)}_{\rm
P_{b}}\right),\quad i=0,\ldots,N$ (30)

and in matrix form it becomes

$\sum_{j=0}^{N}a^{D}_{ij}u_{j}=f^{D}_{i},\quad i=0,\ldots,N$ (31)

where $A^{D}=\\{a^{D}_{i,j}\\}$ and $F^{D}=\\{f^{D}_{i}\\}$, and

$\begin{array}[]{c}\displaystyle A^{D}={\rm S+S_{T}+}\lambda\,\left({\rm
P_{a}+P_{b}}\right)\\\ \displaystyle F^{D}={\rm M}f+\left({\rm
D_{a}}+\lambda{\rm P_{a}}\right)u_{a}+\left({\rm D_{b}}+\lambda{\rm
P_{b}}\right)u_{b}\end{array}$ (32)

where the matrices ${\rm
A^{D},M,D_{a},D_{b},P_{a},P_{b}}\in\mathbb{R}^{(N+1)\times(N+1)}$.

With reference to Fig. 1, we compute the entries of the matrices previously
described to obtain the following values for all non zero entries of the
matrices:

$\begin{array}[]{c}\displaystyle(\varphi_{0}^{\prime},\varphi_{0}^{\prime})_{L^{2}([a,b])}=\frac{\theta_{1}}{h},\quad(\varphi_{0}^{\prime},\varphi_{1}^{\prime})_{L^{2}([a,b])}=-\frac{\theta_{1}}{h},\quad(\varphi_{1}^{\prime},\varphi_{1}^{\prime})_{L^{2}([a,b])}=\frac{1+\theta_{1}}{h},\\\\[11.38109pt]
\displaystyle(\varphi_{N-1}^{\prime},\varphi_{N-1}^{\prime})_{L^{2}([a,b])}=\frac{1+\theta_{2}}{h},\quad(\varphi_{N-1}^{\prime},\varphi_{N}^{\prime})_{L^{2}([a,b])}=-\frac{\theta_{2}}{h},\quad(\varphi_{N}^{\prime},\varphi_{N}^{\prime})_{L^{2}([a,b])}=\frac{\theta_{2}}{h},\\\\[11.38109pt]
\displaystyle(\varphi_{i-1}^{\prime},\varphi_{i}^{\prime})_{L^{2}([a,b])}=-\frac{1}{h},\>i=1,\ldots,N\quad(\varphi_{i}^{\prime},\varphi_{i}^{\prime})_{L^{2}([a,b])}=\frac{2}{h},i=2,\ldots,N-2,\\\\[11.38109pt]
\displaystyle\varphi_{0}(a)=\theta_{1},\quad\varphi_{1}(a)=1-\theta_{1},\quad\varphi_{N}(b)=\theta_{2},\quad\varphi_{N-1}(b)=1-\theta_{2},\\\\[11.38109pt]
\displaystyle\varphi^{\prime}_{0}(a)=-1/h,\quad\varphi^{\prime}_{1}(a)=1/h,\quad\varphi^{\prime}_{N}(b)=1/h,\quad\varphi^{\prime}_{N-1}(b)=-1/h.\end{array}$
(33)

The rest of the source terms can be computed as follows. The term related to
the mass matrix is

$\begin{array}[]{c}\displaystyle(f,\varphi_{i})_{L^{2}([a,b])}\approx\sum_{j=0}^{N}f_{j}(\varphi_{j},\varphi_{i})_{L^{2}([a,b])}=\frac{h}{6}\left(f_{i-1}+4f_{i}+f_{i+1}\right),\quad
i=2,\ldots,n-2.\\\
\displaystyle(f,\varphi_{0})_{L^{2}([a,b])}\approx\sum_{j=0}^{N}f_{j}(\varphi_{j},\varphi_{0})_{L^{2}([a,b])}=h\left(\frac{\theta_{1}^{3}}{3}f_{0}+\left(\frac{\theta_{1}^{2}}{2}-\frac{\theta_{1}^{3}}{3}\right)f_{1}\right),\\\\[14.22636pt]
\displaystyle(f,\varphi_{1})_{L^{2}([a,b])}\approx\sum_{j=0}^{N}f_{j}(\varphi_{j},\varphi_{1})_{L^{2}([a,b])}=h\left(\left(\frac{\theta_{1}^{2}}{2}-\frac{\theta_{1}^{3}}{3}\right)f_{0}+\frac{1+\theta_{1}^{3}}{3}f_{1}\right),\\\\[14.22636pt]
\displaystyle(f,\varphi_{N-1})_{L^{2}([a,b])}\approx\sum_{j=0}^{N}f_{j}(\varphi_{j},\varphi_{N-1})_{L^{2}([a,b])}=h\left(\left(\frac{1+\theta_{2}^{3}}{3}-\theta_{2}^{2}+\theta_{2}\right)f_{N-1}+\left(\frac{\theta_{2}^{2}}{2}-\frac{\theta_{2}^{3}}{3}\right)f_{N}\right)\\\\[14.22636pt]
\displaystyle(f,\varphi_{N})_{L^{2}([a,b])}\approx\sum_{j=0}^{N}f_{j}(\varphi_{j},\varphi_{N-1})_{L^{2}([a,b])}=h\left(\left(\frac{\theta_{2}^{2}}{2}-\frac{\theta_{2}^{3}}{3}\right)f_{N-1}+\frac{\theta_{2}^{3}}{3}f_{N}\right).\end{array}$
(34)

Regarding the penalisation term

$\begin{array}[]{c}\displaystyle
u_{a}\varphi_{0}(a)\approx\sum_{j=0}^{N}{u_{a}}_{j}\varphi_{j}(a)\varphi_{0}(a)=\theta_{1}^{2}{u_{a}}_{0}+\theta_{1}(1-\theta_{1}){u_{a}}_{1}\\\
\displaystyle
u_{a}\varphi_{1}(a)\approx\sum_{j=0}^{N}{u_{a}}_{j}\varphi_{j}(a)\varphi_{1}(a)=\theta_{1}(1-\theta_{1}){u_{a}}_{0}+(1-\theta_{1}^{2}){u_{a}}_{1},\end{array}$
(35)

analogously, we find the coefficients for the penalisation term in $x=b$.

We now turn our attention to the mixed Dirichlet-Neumann problem (17). We
consider the penalised variational formulation in (25) and look for a solution
in $V_{h}$. Replacing $u$ by $u_{h}$ and $v$ by $\varphi_{j}$ in (27), we
obtain the system

$\displaystyle\int_{a}^{b}u_{h}^{\prime}(x)\varphi_{i}^{\prime}(x)\,dx+u_{h}^{\prime}(a)\varphi_{i}(a)+u_{h}(a)\varphi_{i}^{\prime}(a)+\lambda
u_{h}(a)\varphi_{i}(a)$ (36)
$\displaystyle=\int_{a}^{b}f_{h}(x)\varphi_{i}(x)\,dx+{u_{a}}_{h}\varphi_{i}^{\prime}(a)+{g_{b}}_{h}\varphi_{i}(b)+\lambda{u_{a}}_{h}\varphi_{i}(a),\quad\forall
i=1,\cdots,N$

and if we use the expression in (27) for the function $u_{h}$, it becomes

$\displaystyle\sum_{j=0}^{N}u_{j}\underbrace{(\varphi_{j}^{\prime},\varphi_{i}^{\prime})_{L^{2}([a,b])}}_{\rm
S}+\sum_{j=0}^{N}u_{j}\underbrace{\left[\varphi_{j}^{\prime}(a)\varphi_{i}(a)+\varphi_{j}(a)\varphi_{i}^{\prime}(a)\right]}_{\rm
S_{T}}+\lambda\sum_{j=0}^{N}u_{j}\underbrace{\varphi_{j}(a)\varphi_{i}(a)}_{\rm
P_{a}}$
$\displaystyle=\sum_{j=0}^{N}f_{j}\underbrace{(\varphi_{j},\varphi_{i})_{L^{2}([a,b])}}_{\rm
M}\,dx+\sum_{j=0}^{N}{u_{a}}_{j}\underbrace{\varphi_{j}\varphi_{i}^{\prime}(a)}_{\rm
D_{a}}+\sum_{j=0}^{N}{g_{b}}_{j}\underbrace{\varphi_{j}(b)\varphi_{i}(b)}_{\rm
N_{b}}$
$\displaystyle+\lambda\sum_{j=0}^{N}{u_{a}}_{j}\underbrace{\varphi_{j}\varphi_{i}(a)}_{\rm
P_{a}},\quad i=0,\ldots,N$ (37)

which leads to the following system of equations for $u_{i}$:

$\sum_{j=0}^{N}a^{M}_{ij}u_{j}=f^{M}_{i},\quad i=0,\ldots,N$ (38)

where $A^{M}=\\{a^{M}_{i,j}\\}$ and $F^{M}=\\{f^{M}_{i}\\}$, and

$\begin{array}[]{c}\displaystyle A^{M}={\rm S+S_{T}+}\lambda\,{\rm P_{a}}\\\
\displaystyle F^{M}={\rm M}f+\left({\rm D_{a}}+\lambda{\rm
P_{a}}\right)u_{a}+{\rm N_{b}}\,g_{b}\end{array}$ (39)

The non zero values of the various terms are computed as in Eq. (33).

##  4 The two-dimensional problem

Let us consider now the Poisson equation with Dirichlet boundary conditions in
$\Omega\subseteq R\subset\mathbb{R}^{2}$, where $R$ is a rectangular region:

$-\Delta u=f\quad{\rm in}\>\Omega,\quad u=g_{D}\>{\rm on}\>\partial\Omega$
(40)

Multiplying by a test function $v\in H^{1}(\Omega)$, integrating in $\Omega$
and then integrating by parts, one obtains:

$\int_{\Omega}\nabla u\cdot\nabla
v\,d\Omega=\int_{\Omega}fv\,d\Omega+\int_{\partial\Omega}v\frac{\partial
u}{\partial n}\,ds.$ (41)

Following the same strategy that we applied in the one dimensional case, we
penalise the distance from the prescribed value of the Dirichlet boundary
condition by minimizing the following functional

$J_{D\lambda}[u]=\int_{\Omega}\left[\frac{1}{2}(\nabla
u)^{2}-fu\right]\,d\Omega+\frac{\lambda}{2}\int_{\partial\Omega}(u-g_{D})^{2}\,ds$
(42)

By setting to zero the first variation of the functional leads to the
following penalised variational formulation of the problem:

$\int_{\Omega}\nabla u\cdot\nabla
v\,d\Omega+\int_{\partial\Omega}\left[\lambda\cdot(u-g_{D})-\frac{\partial
u}{\partial n}\right]v\,ds=\int_{\Omega}fv\,d\Omega$ (43)

The discretization of the problem is obtained by discretizing the rectangular
region $R\subset\mathbb{R}^{2}$, containing the domain $\Omega$, by a regular
rectangular grid (see Figure 3).

\begin{overpic}[abs,width=158.99156pt,unit=1mm,scale={.25}]{Figures/domain_bc.pdf}
\put(10.0,43.0){$R$} \put(34.5,44.0){$\Gamma_{D}$}
\put(38.0,8.5){$\Gamma_{N}$} \put(28.0,30.0){$\Omega$} \end{overpic} Figure 2:
Two dimensional problem domain $\Omega\subset R$: $\Gamma_{D}$ and
$\Gamma_{N}$ are portions of $\Gamma=\partial\Omega$, with Dirichlet and
Neumann boundary conditions, respectively. $R$ is the rectangular region where
we define the square grid.

Following the approach showed in [45, 39, 40, 43], the domain $\Omega$ is
implicitly defined by a level set function $\phi(x,y)$ that is negative inside
$\Omega$, positive outside and zero on the boundary $\Gamma$ :

$\displaystyle\Omega=\\{(x,y):\phi(x,y)<0\\},\qquad\Gamma=\\{(x,y):\phi(x,y)=0\\}.$
(44)

The unit normal vector $n$ in (63) can be computed as
$n=\frac{\nabla\phi}{|\nabla\phi|}$. For example, for a circular domain
centered in $(x_{c},y_{c})$, the most convenient level-set function in terms
of numerical stability is the signed distance function between $(x,y)$ and
$\Gamma$, i.e. $\phi=r-\sqrt{(x-x_{c})^{2}+(y-y_{c})^{2}}$, where $r$ is the
radius of the circle.

The set of grid points will be denoted by $\mathcal{N}$, with
$\\#\mathcal{N}=(1+N)^{2}$, the active nodes (i.e. internal $\mathcal{I}$ or
ghost $\mathcal{G}$) by
$\mathcal{I}\cup\mathcal{G}={\mathcal{A}}\subset\mathcal{N}$, while the set of
inactive points by $\mathcal{O}\subset\mathcal{N}$, with
$\mathcal{O}\cup\mathcal{A}=\mathcal{N}$ and
$\mathcal{O}\cap\mathcal{A}=\emptyset$ and the set of cells by $\mathcal{C}$,
with $\\#\mathcal{C}=N^{2}$.

Here we define the set of ghost points $\mathcal{G}$, which are points that
belong to $R\setminus\Omega$, with at least an internal point as neighbor, and
are formally defines as follows

$(x_{i},y_{i})\in\mathcal{G}\iff(x_{i},y_{i})\in R\setminus\Omega\text{ and
}\\{(x_{i}\pm h,y_{i}),(x_{i},y_{i}\pm h)\\}\cap\mathcal{I}\neq\emptyset.$
(45)

The discrete space $V_{h}$ is given by the piecewise bilinear functions which
are continuous in $R$. As a basis we choose the following functions:

$\varphi_{i}(x,y)=\max\left\\{\left(1-\frac{|x-x_{i}|}{\Delta
x}\right)\left(1-\frac{|y-y_{i}|}{\Delta y}\right),0\right\\},$ (46)

here we denote by $i=(i_{1},i_{2})$ an index that identifies a node on the
grid.

The generic element $u_{h}\in V_{h}$ will have the following representation:

$u_{h}(x,y)=\sum_{i\in{\mathcal{A}}}u_{i}\varphi_{i}(x,y).$ (47)

The discretization of the problem (13) is obtained by replacing $u$ and $v$ by
$u_{h}$ and $v_{h}$, both in $V_{h}$, and evaluating the integral instead of
over $\Omega$ over its polygonal approximation $\Omega_{h}$. Since we are
approximating $\Omega$ with $\Omega_{h}$, we are also approximating
$\Gamma_{\bullet}$ by its polygonal approximation $\Gamma_{\bullet,h}$. This
leads to the following discrete variational problem

$\int_{\Omega_{h}}\nabla u_{h}\cdot\nabla
v_{h}\,d{\Omega_{h}}+\int_{\Gamma_{h}}\left[\lambda\cdot(u_{h}-g_{D})-\frac{\partial
u_{h}}{\partial
n}\right]v_{h}\,ds=\int_{\Omega_{h}}fv_{h}\,d{\Omega_{h}},\quad\forall
v_{h}\in V_{h}$ (48)

When applied to the subspace $V_{h}$, and using the standard basis, the above
equations become

$\displaystyle\sum_{j\in\mathcal{N}}u_{j}(\nabla\varphi_{i},\nabla\varphi_{j})_{L^{2}(\Omega_{h})}-\sum_{j\in\mathcal{N}}u_{j}\left(\frac{\partial\varphi_{i}}{\partial
n},\varphi_{j}\right)_{L^{2}(\Gamma_{D,h})}+\lambda\sum_{j\in\mathcal{N}}u_{j}(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{D,h})}$
$\displaystyle=$
$\displaystyle\sum_{j\in\mathcal{N}}f_{j}(\varphi_{j},\varphi_{i})_{L^{2}(\Omega_{h})}+\lambda\sum_{j\in\mathcal{N}}{g_{D}}_{j}(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{D,h})}.$
(49)

where the quantities $f_{i}$, $g^{D}_{i}$ denote the node values of source and
of the Dirichlet boundary functions.

Here we rewrite (48), applying the Nitsche technique, to make the equation
symmetric

$\displaystyle\int_{\Omega_{h}}\nabla u_{h}\cdot\nabla
v_{h}\,d{\Omega_{h}}+\int_{\Gamma_{h}}\lambda\cdot(u_{h}-g_{D})v_{h}\,ds-\int_{\Gamma_{h}}\frac{\partial
u_{h}}{\partial n}v_{h}\,ds-\int_{\Gamma_{h}}u_{h}\frac{\partial
v_{h}}{\partial n}\,ds$
$\displaystyle=\int_{\Omega_{h}}fv_{h}\,d{\Omega_{h}}-\int_{\Gamma_{h}}g_{D}\frac{\partial
v_{h}}{\partial n}\,ds,\quad\forall v_{h}\in V_{h}$ (50)

as follows

$\displaystyle\sum_{j\in\mathcal{N}}u_{j}\left(\underbrace{(\nabla\varphi_{j},\nabla\varphi_{i})_{L^{2}(\Omega_{h})}}_{\rm
S}-\underbrace{\left(\left(\frac{\partial\varphi_{j}}{\partial
n},\varphi_{i}\right)_{L^{2}(\Gamma_{D,h})}+\left(\varphi_{j},\frac{\partial\varphi_{i}}{\partial
n}\right)_{L^{2}(\Gamma_{D,h})}\right)}_{\rm
S_{T}}+\lambda\underbrace{(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{D,h})}}_{\rm
P_{a}}\right)$ $\displaystyle=$
$\displaystyle\sum_{j\in\mathcal{N}}f_{j}\underbrace{(\varphi_{j},\varphi_{i})_{L^{2}(\Omega_{h})}}_{\rm
M}+\sum_{j\in\mathcal{N}}{g_{D}}_{j}\left(-\underbrace{\left(\varphi_{j},\frac{\partial\varphi_{i}}{\partial
n}\right)_{L^{2}(\Gamma_{D,h})}}_{\rm
D_{a}}+\underbrace{\lambda(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{D,h})}}_{\rm
P_{a}}\right)$ (51)

Making use of representation (47) and using the functions of the basis as test
functions we obtain the following system of equations for $u_{i}$,
$i\in{\mathcal{A}}$:

$\sum_{j\in{\mathcal{A}}}a^{D}_{ij}u_{j}=f^{D}_{i},\quad i\in{\mathcal{A}}$
(52)

where

$\begin{array}[]{c}\displaystyle
a^{D}_{ij}=(\nabla\varphi_{i},\nabla\varphi_{j})_{L^{2}(\Omega_{h})}+\lambda(\varphi_{i},\varphi_{j})_{L^{2}(\Gamma_{D,h})}-(n\cdot\nabla\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{D,h})}\\\
-(\varphi_{j},n\cdot\nabla\varphi_{i})_{L^{2}(\Gamma_{D,h})},\\\\[11.38109pt]
\displaystyle
f^{D}_{i}=\lambda(g_{D},\varphi_{i})_{L^{2}(\Gamma_{D,h})}-(g_{D},n\cdot\nabla\varphi_{i})_{L^{2}(\Gamma_{D,h})}+(f,\varphi_{i})_{L^{2}(\Omega_{h})}.\end{array}$
(53)

As discussed in the one-dimensional case, it is easy to obtain a variational
formulation for

$-\Delta u=f\quad{\rm in}\>\Omega,\qquad\frac{\partial u}{\partial
n}=g_{N}\>{\rm on}\>\Gamma_{N},$ (54)

when $\Gamma_{N}$ is the all perimeter of $\Omega$, i.e
$\Gamma_{N}=\partial\Omega$. It is sufficient to multiply the previous
equation by a test function $v\in H^{1}(\Omega)$, integrate over $\Omega$ and
then integrating by parts to obtain

$\displaystyle\int_{\Omega}\nabla u\cdot\nabla
v\,d\Omega=\int_{\Omega}fv\,d\Omega+\int_{\Gamma_{N}}g_{N}v\,ds$ (55)

Once again we need to be careful that a compatibility condition has to be
verified by $f\in L^{2}(\Omega)$ and we need to choose a unique solution in
the kernel of the gradient, imposing a null-average condition, i.e.

$\int_{\Omega}f\,d\Omega=\int_{\Gamma_{N}}g_{N}\,ds,\qquad\int_{\Omega}u\,d\Omega=0.$
(56)

Here we discretize (55) using the standard basis of the discrete space
$V_{h}$, that becomes

$\displaystyle\sum_{j\in\mathcal{N}}u_{j}(\nabla\varphi_{i},\nabla\varphi_{j})_{L^{2}(\Omega_{h})}=\sum_{j\in\mathcal{N}}f_{j}(\varphi_{j},\varphi_{i})_{L^{2}(\Omega_{h})}+\sum_{j\in\mathcal{N}}{g_{N}}_{j}(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{N,h})}$
(57)

where the quantities $f_{i}$, $g^{N}_{i}$ denote the node values of source and
of the Neumann boundary functions, respectively. Here we rewrite (57), as
follows

$\displaystyle\sum_{j\in\mathcal{N}}u_{j}\underbrace{(\nabla\varphi_{j},\nabla\varphi_{i})_{L^{2}(\Omega_{h})}}_{\rm
S}=\sum_{j\in\mathcal{N}}f_{j}\underbrace{(\varphi_{j},\varphi_{i})_{L^{2}(\Omega_{h})}}_{\rm
M}+\sum_{j\in\mathcal{N}}{g_{N}}_{j}\underbrace{(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{N,h})}}_{\rm
N_{b}}$

This results is the discrete system

$\sum_{j\in{\mathcal{A}}}a^{N}_{ij}u_{j}=f^{N}_{i},\quad i\in{\mathcal{A}}$
(58)

where $A^{N}=\\{a^{N}_{ij}\\}$ and $F^{N}=\\{f_{i}\\}$, and

$\displaystyle A^{N}={\rm S},$ (59) $\displaystyle F^{N}={\rm M}f+{\rm
N_{b}}\,g_{N},$ (60)

where the entries of the matrices are

$\displaystyle\displaystyle a^{N}_{ij}$
$\displaystyle=(\nabla\varphi_{i},\nabla\varphi_{j})_{L^{2}(\Omega_{h})},$
(61) $\displaystyle\displaystyle f^{N}_{i}$
$\displaystyle=(f,\varphi_{i})_{L^{2}(\Omega_{h})}+(g_{N},\varphi_{i})_{L^{2}(\Gamma_{N,h})}.$

Let us now consider the following mixed problem:

$-u^{\prime\prime}(x)=f(x),\quad x\in[a,b]\qquad u(a)=0,\quad
u^{\prime}(b)=g_{b}$ (62)

We now consider the mixed problem in two dimensions:

$-\Delta u=f\quad{\rm in}\>\Omega,\quad u=g_{D}\>{\rm
on}\>\Gamma_{D},\>\frac{\partial u}{\partial n}=g_{N}\>{\rm on}\>\Gamma_{N}$
(63)

where the boundary $\partial\Omega$ is partitioned into a Dirichlet boundary
$\Gamma_{D}$ and a Neumann boundary $\Gamma_{N}$,
$\Gamma_{D}\cup\Gamma_{N}=\partial\Omega$,
$\Gamma_{D}\cap\Gamma_{N}=\emptyset$.

The penalised method is obtained by minimizing the following functional, which
differs from the original one because of the addition of a penalisation term
on the Dirichlet boundary $\Gamma_{D}$:

$J_{D\lambda}[u]=\int_{\Omega}\left[\frac{1}{2}(\nabla
u)^{2}-fu\right]\,d\Omega+\frac{\lambda}{2}\int_{\Gamma_{D}}(u-g_{D})^{2}\,ds$
(64)

Setting to zero the first variation of the functional leads to the following
penalised variational formulation of the problem:

$\displaystyle\int_{\Omega}\nabla u\cdot\nabla
v\,d\Omega+\int_{\Gamma_{D}}\Big{[}\lambda\cdot(u-g_{D})$
$\displaystyle-\frac{\partial u}{\partial n}\Big{]}v\,ds$ $\displaystyle=$
$\displaystyle\int_{\Omega}fv\,d\Omega+\int_{\Gamma_{N}}g_{N}v\,ds$ (65)

where we made use of Neumann boundary conditions by replacing
$\displaystyle\frac{\partial u}{\partial n}$ with $g_{N}$ on $\Gamma_{N}$.
Also in two space dimensions it is possible to make symmetric the problem by
adopting the Nitsche technique, adding and subtracting the same quantity, as
follows

$\displaystyle\int_{\Omega}\nabla u\cdot\nabla v\,d\Omega+$
$\displaystyle\int_{\Gamma_{D}}\Big{[}\lambda\cdot(u-g_{D})-\frac{\partial
u}{\partial n}\Big{]}v\,ds-\int_{\Gamma_{D}}u\frac{\partial v}{\partial
n}\,ds$ $\displaystyle=$
$\displaystyle\int_{\Omega}fv\,d\Omega+\int_{\Gamma_{N}}g_{N}v\,ds-\int_{\Gamma_{D}}g_{D}\frac{\partial
v}{\partial n}\,ds.$ (66)

Notice that for purely Neumann problems there is no penalisation term, since
$\Gamma_{N}=\Gamma$ and $\Gamma_{D}=\emptyset$. Of course for such Neumann
problem we assume the standard compatibility condition

$\int_{\Omega}f\,d\Omega+\int_{\Gamma}g_{N}\,d\Gamma=0.$ (67)

The discretization of the problem is obtained in the usual way, by replacing
the domain $\Omega$ with its polygonal approximation $\Omega_{h}$, whose
boundary is $\Gamma_{h}$.

We now proceed to consider the variational formulation (66) over the discrete
space $V_{h}$, and using the standard basis, the above equations become

$\displaystyle\sum_{j\in\mathcal{N}}u_{j}(\nabla\varphi_{i},\nabla\varphi_{j})_{L^{2}(\Omega_{h})}-\sum_{j\in\mathcal{N}}u_{j}\left(\frac{\partial\varphi_{i}}{\partial
n},\varphi_{j}\right)_{L^{2}(\Gamma_{D,h})}-\sum_{j\in\mathcal{N}}u_{j}\left(\varphi_{j},\frac{\partial\varphi_{i}}{\partial
n}\right)_{L^{2}(\Gamma_{D,h})}$ $\displaystyle+$
$\displaystyle\lambda\sum_{j\in\mathcal{N}}u_{j}(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{D,h})}=\sum_{j\in\mathcal{N}}f_{j}(\varphi_{j},\varphi_{i})_{L^{2}(\Omega_{h})}+\sum_{j\in\mathcal{N}}{g_{N}}_{j}(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{N,h})}$
(68) $\displaystyle-$
$\displaystyle\sum_{j\in\mathcal{N}}{g_{D}}_{j}\left(\varphi_{j},\frac{\partial\varphi_{i}}{\partial
n}\right)_{L^{2}(\Gamma_{D,h})}+\lambda\sum_{j\in\mathcal{N}}{g_{D}}_{j}(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{D,h})}$

where the quantities $f_{i}$, $g^{D}_{i}$, $g^{N}_{i}$ denote the node values
of source and of the Dirichlet and Neumann boundary functions, respectively.
Here we rewrite (4), as follows

$\displaystyle\sum_{j\in\mathcal{N}}u_{j}\left(\underbrace{(\nabla\varphi_{j},\nabla\varphi_{i})_{L^{2}(\Omega_{h})}}_{\rm
S}-\underbrace{\left(\left(\frac{\partial\varphi_{j}}{\partial
n},\varphi_{i}\right)_{L^{2}(\Gamma_{D,h})}+\left(\varphi_{j},\frac{\partial\varphi_{i}}{\partial
n}\right)_{L^{2}(\Gamma_{D,h})}\right)}_{\rm S_{T}}\right.$ $\displaystyle+$
$\displaystyle\left.\lambda\underbrace{(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{D,h})}}_{\rm
P_{a}}\right)=\sum_{j\in\mathcal{N}}f_{j}\underbrace{(\varphi_{j},\varphi_{i})_{L^{2}(\Omega_{h})}}_{\rm
M}+\sum_{j\in\mathcal{N}}{g_{N}}_{j}\underbrace{(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{N,h})}}_{\rm
N_{b}}$ (69)
$\displaystyle+\sum_{j\in\mathcal{N}}{g_{D}}_{j}\left(-\underbrace{\left(\varphi_{j},\frac{\partial\varphi_{i}}{\partial
n}\right)_{L^{2}(\Gamma_{D,h})}}_{\rm
D_{a}}\\!\\!\\!\\!\\!\\!\\!+\underbrace{\lambda(\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{D,h})}}_{\rm
P_{a}}\right)$ (70)

This results is the discrete system

$\sum_{j\in{\mathcal{A}}}a^{M}_{ij}u_{j}=f^{M}_{i},\quad i\in{\mathcal{A}}$
(71)

where $A^{M}=\\{a^{M}_{ij}\\}$ and $F^{M}=\\{f_{i}\\}$, and

$\displaystyle A^{M}={\rm S-S_{T}}+\lambda{\rm P_{a}},$ (72) $\displaystyle
F^{M}={\rm M}f+(\lambda{\rm P_{a}-D_{a}})g_{D}+{\rm N_{b}}\,g_{N},$

where the entries of the matrices are

$\begin{array}[]{c}\displaystyle
a^{M}_{ij}=(\nabla\varphi_{i},\nabla\varphi_{j})_{L^{2}(\Omega_{h})}+\lambda(\varphi_{i},\varphi_{j})_{L^{2}(\Gamma_{D,h})}-(n\cdot\nabla\varphi_{j},\varphi_{i})_{L^{2}(\Gamma_{D,h})}\\\
-(\varphi_{j},n\cdot\nabla\varphi_{i})_{L^{2}(\Gamma_{D,h})},\\\\[11.38109pt]
\displaystyle
f^{M}_{i}=\lambda(g_{D},\varphi_{i})_{L^{2}(\Gamma_{D,h})}-(g_{D},n\cdot\nabla\varphi_{i})_{L^{2}(\Gamma_{D,h})}+(g_{N},\varphi_{i})_{L^{2}(\Gamma_{N,h})}\\\
+(f,\varphi_{i})_{L^{2}(\Omega_{h})}.\end{array}$ (73)

##  5 Computational aspects and snapping back to grid

As previously discussed we are interested in dealing with the discrete
variational formulation (50) and (4). For this reason, in this section, we
will detail the implementation aspects regarding how to compute the entries of
the linear system corresponding to (50) and (4).

In particular, we will adopt a reference element type of approach that is
standard in the finite element community. Hence, we consider a structured
tessellation $\mathcal{T}_{h}$ (see Fig. 4 (a)) of the computational domain
$R\subset\mathbb{R}^{2}$ and map each element of the tessellation to reference
unit square $\widehat{K}$ where (50) and (4) are evaluated, [11, 18].

In particular, we will denote $F_{K}:K\to\widehat{K}$ the map from the
physical element $K\in\mathcal{T}_{h}$ to the reference element $\widehat{K}$.

\begin{overpic}[abs,width=278.23827pt,unit=1mm,scale={.25}]{Figures/Domain_mesh.pdf}
\put(-7.0,58.0){(a)} \put(-2.0,44.0){$\nearrow$}
\put(-2.0,25.0){$\xrightarrow{\hskip 56.9055pt}$} \put(-15.0,27.0){$\phi<0$}
\put(-15.0,40.0){$\phi>0$} \put(48.0,58.0){(b)} \put(21.5,30.8){$\Omega$}
\put(11.0,56.0){$R$} \put(65.0,56.0){$R$} \put(80.0,35.0){$\Omega$}
\end{overpic} Figure 3: (a): Square regular mesh of the region R and of the
domain $\Omega$. In red we have those cells that intersect the boundary
$\Gamma$. (b): Representation of the two points A, B defined in Algorithm 1,
making use of the points $\\{P_{0},\ldots,P_{4}\\}$, where $P_{0}$ and $P_{1}$
are ghost points.

\begin{overpic}[abs,width=317.9892pt,unit=1mm,scale={.25}]{Figures/polygon_shapes_border.pdf}
\put(-5.0,48.0){(a)} \put(27.0,47.0){convex} \put(48.0,32.0){non--convex}
\put(28.0,20.0){$\Omega$} \end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/index_local_global.pdf}
\put(-5.0,40.0){(b)} \put(18.5,5.5){$\searrow$} \put(25.0,28.0){$\Omega$}
\put(52.0,16.0){$\widehat{K}$} \put(44.0,13.0){$\widehat{l}_{0}$}
\put(53.0,25.0){$\widehat{l}_{1}$} \put(53.0,8.5){$\widehat{l}_{4}$}
\put(65.5,17.0){$\widehat{l}_{2}$} \put(60.0,6.0){$\searrow$}
\put(65.0,2.0){$\widehat{l}_{3}$} \end{overpic}

Figure 4: (a): Possible shapes of those polygons that intersect the boundary
$\Gamma$. These elements can assume different shapes, from triangles to
quadrilateral and pentagons. In red we show the set of cells $\mathcal{C}$
after the snapping back to grid. (b): Representation of the indices
$\\{k_{0},\ldots,k_{m-1}\\}$, the edges $\\{l_{0},\ldots,l_{m-1}\\}$ of a
generic element $K$, and the intersections with $\Gamma$, $\bf A$ and $\bf B$;
on the right, a zoom-in that represents the reference element $\widehat{K}$,
with edges $\\{\widehat{l}_{0},\ldots,\widehat{l}_{m-1}\\}$, and vertices
$\\{P_{1},\ldots,P_{m+1}\\}$.

If $i$ is an inactive point, i.e. $i\not\in\mathcal{A}$, we shall set
$f^{\bullet}_{i}=0$, and the components of matrix $A^{\bullet}$ corresponding
to such a point are $a^{\bullet}_{ij}=a^{\bullet}_{ji}=\delta_{ij}$. In this
way matrix $A^{\bullet}$ will be non singular and the solution will be zero on
the inactive points. We adopted the notation
$f^{\bullet},a^{\bullet}_{ij},A^{\bullet}$ to denote either
$f^{D},a^{D}_{ij},A^{D}$, $f^{N},a^{N}_{ij},A^{N}$ or $f^{M},a^{M}_{ij},A^{M}$
as in (39),(61) and (72).

We consider now the various integrals appearing in (53),(61) and (73). Let us
start from the integral $\left(\varphi_{i},\varphi_{j}\right)_{L^{2}(K)}$.
Since $\varphi_{{j}}({x}_{i},y_{i})=\delta_{{i}{j}},\,i,j\in\mathcal{N}$, it
is possible to reduce the summations to the vertices of the $K-$th cell (see
Fig. 4 (b)), yielding

$\displaystyle\left(\varphi_{i},\varphi_{j}\right)_{L^{2}(K)}=\sum_{\eta,\mu=0}^{m-1}\left(\varphi_{k_{\eta}},\varphi_{k_{\mu}}\right)_{L^{2}(K)}.$
(74)

where we denote by capital letters the cell indices, i.e. $K\in\mathcal{C}$,
and by small letters the vertices of the cells.

Algorithm 1 Computation of the intersection of the boundary with the grid

$k_{4}=k_{0}$

for i = 0:3 do

if $\phi(k_{i})\phi(k_{i+1})<0$ then

$\theta=\phi(k_{i})/(\phi(k_{i})-\phi(k_{i+1}))$

$P=\theta k_{i+1}+(1-\theta)k_{i}$

if $\phi(k_{i})<0$ then

${\bf A}:=P$

else

${\bf B}:=P$

end if

end if

end for

Two fundamental novelties of the scheme we propose are the following: first,
the unknowns are always located on the original square grid, more precisely,
on interior or ghost nodes. Second, the computation of the integrals that
define all matrix elements, characterizing systems (50) and (4), are computed
by exact quadrature formulas defined on the boundary of each cell, as
explained below.

We observe that the product of two elements in $V_{h}$ appearing in (74), is
an element of $\mathbb{Q}_{2}(K)$, i.e. the set of bi-quadratic polynomials in
$K$.

Considering a generic function $f$ defined in $\Omega$, if we chose $F=\int
fdx$, it is $\nabla\cdot\mathbf{F}=f$, where $\mathbf{F}=(F,0)^{\top}$. Thus
we have

$\displaystyle\int_{K}fdK=\int_{K}\nabla\cdot\mathbf{F}\,dK.$ (75)

We can apply the Gauss theorem to (75) to obtain

$\displaystyle\int_{K}\nabla\cdot\mathbf{F}\,dK=\int_{\partial
K}\mathbf{F}\cdot\mathbf{n}\,dl.$ (76)

We now consider $\mathcal{P}$, the generic polygon with $m$ edges
$l_{i},\,i=0,\ldots,m-1$ (see Fig. 4 (b) and Fig. 5) and express (76) as

$\displaystyle\int_{\mathcal{P}}\nabla\cdot\mathbf{F}\,d\mathcal{P}=\int_{\partial\mathcal{P}}\mathbf{F}\cdot\mathbf{n}\,dl=\sum_{i=0}^{m-1}\int_{l_{i}}\mathbf{F}\cdot\mathbf{n}\,dl=\sum_{i=0}^{m-1}\int_{l_{i}}F{n}_{x}\,dl=\sum_{i=0}^{m-1}\int_{l_{i}}Fdy.$
(77)

To evaluate the integral over the generic edge $l_{i}$ we choose the three-
point Gauss-Legendre quadrature rule (which is exact for polynomials in
$\mathbb{P}_{5}(\mathbb{R})$), thus we write

$\int_{l_{i}}Fdy=\sum_{s=1}^{3}w_{s}F(\widehat{x}_{i,s},\widehat{y}_{i,s})(y_{P_{i+1}}-y_{P_{i}})$
(78)

(see Fig. 4 (b)) where $w_{j}\in\\{5/18,4/9,5/18\\}$,
$\widehat{x}_{j}\in\\{1-\sqrt{3/5},1,1+\sqrt{3/5}\\}/2\subset[0,1]$, with
$j=1,2,3$.

Figure 5: Scheme of the three quadrature points (circles) for each edge
$l_{i},\,i=0,\cdots,4$. The squared points represent the vertices
$P_{i},\,i=0,\cdots,4,$ of the polygon $\mathcal{P}$.

Choosing
$f=\varphi_{k_{l}}\varphi_{k_{m}}\in\mathbb{Q}_{2}(K)\subset\mathbb{P}_{4}(K)$
then $F\in\mathbb{P}_{5}(K)$, and making use of (77), we write

$\left(\varphi_{i},\varphi_{j}\right)_{L^{2}(K)}=\sum_{i=0}^{m-1}\left(\sum_{s=1}^{3}w_{s}\Phi_{i,j}(\widehat{x}_{i,s},\widehat{y}_{i,s})(y_{P_{i+1}}-y_{P_{i}})\right)|l_{i}|.$
(79)

where $\Phi_{i,j}=\int\varphi_{i}\varphi_{j}\,dx$.

Analogously, the integral for the stiff term in (53), (73) and (73) is
computed as follows

$\left(\nabla\varphi_{i},\nabla\varphi_{j}\right)_{L^{2}(K)}=\sum_{i=0}^{m-1}\left(\sum_{s=1}^{3}w_{s}\Psi_{i,j}(\widehat{x}_{i,s},\widehat{y}_{i,s})(y_{P_{i+1}}-y_{P_{i}})\right)|l_{i}|.$
(80)

where $\Psi_{i,j}=\int\nabla\varphi_{i}\nabla\varphi_{j}\,dx$.

Stability analysis plays an important role in the analysis of numerical
schemes when solving a partial differential equation in arbitrary domains
[13]. In Fig. 4, we show different types of stability failure, where we modify
the discrete domain $\Omega_{h}$. When the surface of a cell, that intersects
the boundary, is proportional to a power of the space step, i.e. $\propto
h^{\alpha}$, we neglect that cell. These cases are trivial because the
computation of the integrals in (50) and (4) could increase the conditioning
number of the computational matrices $A^{D}$ or $A^{M}$. The snapping back to
grid technique, that we show in Algorithm 2 (see Fig. 4 (a)) is a way to avoid
an ill–conditioned matrix, while maintaining the order of accuracy of the
method.

Various techniques, both physical and mathematical, were employed to discover
the optimal shape of the computational domain that improves the stability of
the numerical schemes [37]. In Algorithm 2 we identify internal cells located
in proximity to the boundary $\Gamma$. Evaluating the level set function
$\phi$ at the vertices of each cell, if the value is less than $h^{\alpha}$,
for a chosen $\alpha$, we disregard the respective cell.

In Fig. 4 (a), we present two distinct scenarios starting from a convex
domain: one, where the remaining computational domain retains its convex
nature, and a second scenario where convexity is lost due to the application
of the snapping back to grid technique. For these two different cases, a
different analysis will be conducted in Section 6.

Algorithm 2 Snapping back to grid

for $k\in\mathcal{N}$ do

if $\phi(k)<0\quad\&\quad|\phi(k)|<h^{\alpha}$ then

$\phi(k):=0$

end if

end for

##  6 A priori convergence analysis

In this section we will focus our attention on developing a priori error
estimates for the scheme described in the paper. The analysis here presented
will follows closely [35]. To lighten the notation we introduce the discrete
bilinear form $a_{h,\alpha}:H^{\frac{3}{2}+\varepsilon}(\Omega)\times
H^{\frac{3}{2}+\varepsilon}(\Omega)\to\mathbb{R}$ corresponding to (4), i.e.

$a_{h,\alpha}(u,v):=(\nabla u,\nabla
v)_{L^{2}(\Omega)}+C_{a}h^{-\alpha}(u,v)_{L^{2}(\Gamma_{D})}-(\partial_{n}u,v)_{L^{2}(\Gamma_{D})}-(u,\partial_{n}v)_{L^{2}(\Gamma_{D})}.$
(81)

Notice that the bilinear form $a_{h}(u,v)$ make sense even when $u$ lives in
$H^{\frac{3}{2}+\varepsilon}_{\Gamma_{D}}$, for this reason we will require
that $(f,g_{D},g_{N})\in L^{2}(\Omega)\times L^{2}(\partial\Omega)\times
H^{-\frac{1}{2}}(\partial\Omega)$ and that $\Omega$ is a polygonal domain. The
request $(f,g_{D},g_{N})\in L^{2}(\partial\Omega)\times
L^{2}(\partial\Omega)\times H^{-\frac{1}{2}}(\partial\Omega)$ and $\Omega$
polygonal are sufficient to ensure that $u\in H^{1}_{\Gamma_{D}}(\Omega)$ such
that

$a(u,v):=(\nabla u,\nabla
v)_{L^{2}(\Omega)}=(f,v)_{L^{2}(\Omega)}+(g_{N},v)_{L^{2}(\Gamma_{N})},$ (82)

lives in $H^{\frac{3}{2}+\varepsilon}_{\Gamma_{D}}$ by elliptic regularity,
[30].

To study the well-posedness of the discrete variational problem, find
$u_{h}\in V_{h}$ such that

$a_{h,\alpha}(u_{h},v_{h}):=(f,v_{h})_{L^{2}(\Omega)}+(g_{N},v_{h})_{L^{2}(\Gamma_{N})}-(g_{D},\partial_{n}v_{h})_{L^{2}(\Gamma_{D})},\qquad
v_{h}\in V_{h}$ (83)

we will consider an auxiliary eigenvalue problem, i.e. find
$(w_{h},\lambda)\in V_{h}\times\mathbb{R}$ such that

$b(w_{h},z_{h}):=(\partial_{n}w_{h},\partial
z_{h})_{L^{2}(\partial\Omega)}=\lambda a(w_{h},z_{h}),\qquad\forall z_{h}\in
V_{h}.$ (84)

We characterize the first eigenvalue problem using the Rayleigh quotient we
have,

$\lambda=\underset{z_{h}\in
V_{h}}{\max}\frac{b(z_{h},z_{h})}{a(z_{h},z_{h})}=\underset{z_{h}\in
V_{h}}{\max}\lVert\partial_{n}z_{h}\rVert_{L^{2}(\partial\Omega)}^{2}\lVert\nabla
z_{h}\rVert^{-2}_{L^{2}(\partial\Omega)}.$ (85)

The above characterization ensures that for any $v_{h}\in V_{h}$ the following
estimate holds

$\lVert\partial_{n}v_{h}\rVert^{2}_{L^{2}(\Omega)}\leq\lambda\lVert\nabla
v_{h}\rVert_{L^{2}(\Omega)}^{2}.$ (86)

We will now bound $\lambda$ in terms of the mesh size, this will be useful to
prove the coercivity of the bilinear form $a_{h}(\cdot,\cdot)$ for
sufficiently large value of $\alpha$. The following lemma is an adaptation of
[27, Lemma 12.8].

######  Lemma 6.1.

Let $V_{h}$ be the space of bilinear functions on each element $K$ of
$\mathcal{T}_{h}$, then

$\lVert\partial_{n}u_{h}\rVert^{2}_{L^{2}(\partial\Omega\cap K)}\leq
C_{2}h^{\gamma-\beta-1}\lVert\nabla u_{h}\rVert_{L^{2}(K)}^{2},\qquad\forall
v_{h}\in V_{h},$ (87)

provided that the following geometric assumption are verified on by
$\mathcal{T}_{h}$:

$\widehat{C}_{\Omega}h^{\beta}\leq\lvert\partial\Omega\cap K\rvert\leq
C_{\Omega}h^{\gamma}.$ (88)

###### Proof.

We will obtain this estimate via a scaling argument, for this reason we
consider the map from the physical element $K$ to the reference element
$\widehat{K}$, i.e. $F_{K}:K\to\widehat{K}$. Notice that since we are working
with a structured mesh the map $F_{K}:K\to\widehat{K}$ is an affine map, i.e.
$F_{K}(\mathbf{x})=\mathbb{A}\mathbf{x}+\mathbf{b}$. We then consider the pull
back of $v_{h}$, i.e.

$\widehat{v}_{h}(\widehat{\mathbf{x}}):=v_{h}\left(F_{K}^{-1}(\widehat{\mathbf{x}})\right):\widehat{K}\to\mathbb{R}.$
(89)

We can then compute via a scaling argument the following bound,

$\displaystyle\lVert\partial_{n}u_{h}\rVert_{L^{2}(\partial\Omega\cap K)}$
$\displaystyle=\lVert\nabla u_{h}\cdot n\rVert_{L^{2}(\partial\Omega\cap
K)}\leq\lVert\mathbf{n}\rVert_{L^{\infty}(\partial\Omega\cap K)}\lVert\nabla
u_{h}\rVert_{L^{2}(\partial\Omega\cap K)}$ (90)
$\displaystyle\leq\lVert\mathbb{A}\rVert^{-1}\left(\frac{\lvert\partial\Omega\cap
K\rvert}{\lvert F_{K}(\partial\Omega\cap
K)\rvert}\right)^{\frac{1}{2}}\lVert\nabla\widehat{u}_{h}\rVert_{L^{2}\left(F_{K}(\partial\Omega\cap
K)\right)}.$ (91)

We now notice that since all norm in finite dimension are equivalent there
exists a constant $C_{1}$ such that

$\lVert\nabla\widehat{u}_{h}\rVert_{L^{2}\left(F_{K}(\partial\Omega\cap
K)\right)}\leq C_{1}\lVert\nabla\widehat{u}_{h}\rVert_{L^{2}(\widehat{K})},$
(92)

where $C_{1}$ may depend on $\beta$ and $\gamma$ but it does not depend on $h$
because we are working on the reference element. Combining (92) with (90) we
have

$\displaystyle\lVert\partial_{n}u_{h}\rVert_{L^{2}(\partial\Omega\cap K)}$
$\displaystyle\leq
C_{1}\lVert\mathbb{A}\rVert^{-1}\left(\frac{\lvert\partial\Omega\cap
K\rvert}{\lvert F_{K}(\partial\Omega\cap
K)\rvert}\right)^{\frac{1}{2}}\lVert\nabla\widehat{u}_{h}\rVert_{L^{2}(\widehat{K})}$
(93) $\displaystyle\leq
C_{1}\lVert\mathbb{A}\rVert^{-1}\lVert\mathbb{A}\rVert\left(\frac{\lvert\partial\Omega\cap
K\rvert}{\lvert F_{K}(\partial\Omega\cap
K)\rvert}\right)^{\frac{1}{2}}\left(\frac{\lvert\widehat{K}\rvert}{\lvert
K\rvert}\right)^{\frac{1}{2}}\lVert\nabla u_{h}\rVert_{L^{2}(K)}$ (94)
$\displaystyle\leq C_{1}\left(\frac{\lvert\widehat{K}\rvert}{\lvert
F_{K}(\partial\Omega\cap
K)\rvert}\right)^{\frac{1}{2}}\left(\frac{\lvert\partial\Omega\cap
K\rvert}{\lvert K\rvert}\right)^{\frac{1}{2}}\lVert\nabla
u_{h}\rVert_{L^{2}(K)}$ (95)

We proceed bounding $\left(\frac{\lvert\widehat{K}\rvert}{\lvert
F_{K}(\partial\Omega\cap
K)\rvert}\right)^{\frac{1}{2}}\left(\frac{\lvert\partial\Omega\cap
K\rvert}{\lvert K\rvert}\right)^{\frac{1}{2}}$ using (88). We begin observing
that $\widehat{K}$ has a unit area and that from (88) we know

$\lvert F_{K}(\partial\Omega\cap
K)\rvert\geq\widehat{C}_{\Omega}h^{\beta-1}\;\Rightarrow\;\lvert
F_{K}(\partial\Omega\cap
K)\rvert^{-1}\leq\widehat{C}_{\Omega}^{-1}h^{1-\beta}.$ (96)

Furthermore $\lvert K\rvert=h^{2}$ and from (88) we know that

$\frac{\lvert\partial\Omega\cap K\rvert}{\lvert K\rvert}\leq
C_{\Omega}h^{\gamma-2}.$ (97)

Combining (96), (97) with (95) we have,

$\displaystyle\lVert\partial_{n}u_{h}\rVert_{L^{2}(\partial\Omega\cap K)}$
$\displaystyle\leq C_{1}\left(\frac{\lvert\widehat{K}\rvert}{\lvert
F_{K}(\partial\Omega\cap
K)\rvert}\right)^{\frac{1}{2}}\left(\frac{\lvert\partial\Omega\cap
K\rvert}{\lvert K\rvert}\right)^{\frac{1}{2}}\lVert\nabla
u_{h}\rVert_{L^{2}(K)}$ (98) $\displaystyle\leq
C_{1}C_{\Omega}^{\frac{1}{2}}\widehat{C}_{\Omega}^{-\frac{1}{2}}h^{\frac{\gamma-\beta-1}{2}}\lVert\nabla
u_{h}\rVert_{L^{2}(K)},$ (99)

lastly we call $C_{2}^{\frac{1}{2}}$ the constant equal to
$C_{1}\widehat{C}_{\Omega}^{-\frac{1}{2}}C_{\Omega}^{\frac{1}{2}}$. ∎

To study the a priori convergence rate of (83) we need to introduce a discrete
norm, i.e.

$\lVert\cdot\rVert_{h,\alpha}^{2}=\lvert\cdot\rvert_{H^{1}(\Omega)}^{2}+h^{-\alpha}C_{a}\lVert\cdot\rVert^{2}_{L^{2}(\Gamma_{D})}+C_{a}^{-1}h^{\alpha}\lVert\partial_{n}\cdot\rVert^{2}_{L^{2}(\Gamma_{D})},$
(100)

which is well posed for any argument in $H^{\frac{3}{2}+\varepsilon}(\Omega)$.

######  Lemma 6.2.

The bilinear form $a_{h,\alpha}:H^{\frac{3}{2}+\varepsilon}(\Omega)\times
H^{\frac{3}{2}+\varepsilon}(\Omega)\to\mathbb{R}$, defined as in (81), is
continuous and coercive, when restricted to $V_{h}$, with respect to
$\lVert\cdot\rVert_{h,\alpha}$, provided we choose $C_{a}=2C_{2}^{-1}$ and
$\alpha\leq\gamma-\beta-1$.

###### Proof.

We begin proving coercivity. We have

$\displaystyle a_{h}(v_{h},v_{h})$ $\displaystyle=(\nabla v_{h},\nabla
v_{h})_{L^{2}(\Omega)}+C_{a}h^{-\alpha}(v_{h},v_{h})_{L^{2}(\Gamma_{D})}-2(\partial_{n}v_{h},v_{h})_{L^{2}(\Gamma_{D})}$
(101) $\displaystyle=\lvert
v_{h}\rvert^{2}_{H^{1}(\Omega)}+C_{a}h^{-\alpha}\lVert
v_{h}\rVert^{2}_{L^{2}(\Gamma_{D})}-2(C_{a}^{\frac{\alpha}{2}}h^{-\frac{\alpha}{2}}\partial_{n}v_{h},C_{a}^{-\frac{\alpha}{2}}h^{\frac{\alpha}{2}}v_{h})_{L^{2}(\Gamma_{D})}$
(102) $\displaystyle\geq\lvert
v_{h}\rvert^{2}_{H^{1}(\Omega)}+(1-\varepsilon)C_{a}h^{-\alpha}\lVert
v_{h}\rVert^{2}_{L^{2}(\Gamma_{D})}-\varepsilon^{-1}C_{a}h^{\alpha}\lVert\partial_{n}v_{h}\rVert^{2}_{L^{2}(\Gamma_{D})},$
(103)

where the last inequality comes from Young’s inequality with
$\varepsilon\in(0,1)$. We add and subtract $\varepsilon\lvert
v_{h}\rvert_{H^{1}(\Omega)}^{2}$, fix $\alpha\leq\gamma-\beta-1$ and use (87)
to control from below $\varepsilon\lvert v_{h}\rvert_{H^{1}(\Omega)}^{2}$,
i.e.

$\displaystyle a_{h,\alpha}(v_{h},v_{h})$ $\displaystyle\geq\lvert
v_{h}\rvert^{2}_{H^{1}(\Omega)}+(1-\varepsilon)C_{a}h^{-\alpha}\lVert
v_{h}\rVert^{2}_{L^{2}(\Gamma_{D})}-\varepsilon^{-1}C_{a}h^{\alpha}\lVert\partial_{n}v_{h}\rVert^{2}_{L^{2}(\Gamma_{D})}$
(104) $\displaystyle\geq(1-\varepsilon)\lvert
v_{h}\rvert^{2}_{H^{1}(\Omega)}+(1-\varepsilon)C_{a}h^{-\alpha}\lVert
v_{h}\rVert^{2}_{L^{2}(\Gamma_{D})}$ (105)
$\displaystyle\qquad\qquad\qquad\qquad+(2\varepsilon-\varepsilon^{-1})C_{a}h^{\alpha}\lVert\partial_{n}v_{h}\rVert^{2}_{L^{2}(\Gamma_{D})}\geq
C_{3}\lVert v_{h}\rVert_{h,\alpha},$ (106)

where we used $C_{a}=2C_{2}^{-1}$, and $C_{3}=(1-\varepsilon^{*})$, with
$\varepsilon^{*}$ being the positive root of
$1-\varepsilon=2\varepsilon-\varepsilon^{-1}$. In particular, notice that
$C_{3}\geq\frac{1}{5}$. We now proceed proving the continuity of the bilinear
form for $u,v\in H^{\frac{3}{2}+\varepsilon}(\Omega)$. We have

$\displaystyle a_{h,\alpha}(u,v)$ $\displaystyle\leq\lvert
u\rvert_{H^{1}(\Omega)}\lvert v\rvert_{H^{1}(\Omega)}$ $\displaystyle+\lVert
h^{-\frac{\alpha}{2}}C_{a}^{\frac{\alpha}{2}}\partial_{n}u\rVert_{L^{2}(\Gamma_{D})}\lVert
h^{\frac{\alpha}{2}}C_{a}^{-\frac{\alpha}{2}}v\rVert_{L^{2}(\Gamma_{D})}$
(107) $\displaystyle+\lVert
h^{-\frac{\alpha}{2}}C_{a}^{\frac{\alpha}{2}}\partial_{n}v\rVert_{L^{2}(\Gamma_{D})}\lVert
h^{\frac{\alpha}{2}}C_{a}^{-\frac{\alpha}{2}}u\rVert_{L^{2}(\Gamma_{D})}$
(108) $\displaystyle+\lVert
h^{-\frac{\alpha}{2}}C_{a}^{\frac{\alpha}{2}}u\rVert_{L^{2}(\Gamma_{D})}\lVert
h^{-\frac{\alpha}{2}}C_{a}^{\frac{\alpha}{2}}v\rVert_{L^{2}(\Gamma_{D})}$
(109) $\displaystyle\leq\lvert u\rvert_{H^{1}(\Omega)}\lvert
v\rvert_{H^{1}(\Omega)}$
$\displaystyle+\left(h^{-\frac{\alpha}{2}}C_{a}^{\frac{\alpha}{2}}\lVert\partial_{n}u\rVert_{L^{2}(\Gamma_{D})}+h^{\frac{\alpha}{2}}C_{a}^{-\frac{\alpha}{2}}\lVert
u\rVert_{L^{2}(\Gamma_{D})}\right)$ (110)
$\displaystyle\cdot\left(h^{-\frac{\alpha}{2}}C_{a}^{\frac{\alpha}{2}}\lVert\partial_{n}v\rVert_{L^{2}(\Gamma_{D})}+h^{\frac{\alpha}{2}}C_{a}^{-\frac{\alpha}{2}}\lVert
v\rVert_{L^{2}(\Gamma_{D})}\right)$ (111) $\displaystyle a_{h,\alpha}(u,v)$
$\displaystyle\leq\left[\lvert
u\rvert^{2}_{H^{1}(\Omega)}+\left(h^{-\frac{\alpha}{2}}C_{a}^{\frac{\alpha}{2}}\lVert\partial_{n}u\rVert_{L^{2}(\Gamma_{D})}+h^{\frac{\alpha}{2}}C_{a}^{-\frac{\alpha}{2}}\lVert
u\rVert_{L^{2}(\Gamma_{D})}\right)^{2}\right]^{\frac{1}{2}}$ (112)
$\displaystyle\cdot\;\;\left[\lvert
v\rvert^{2}_{H^{1}(\Omega)}+\left(h^{-\frac{\alpha}{2}}C_{a}^{\frac{\alpha}{2}}\lVert\partial_{n}v\rVert_{L^{2}(\Gamma_{D})}+h^{\frac{\alpha}{2}}C_{a}^{-\frac{\alpha}{2}}\lVert
v\rVert_{L^{2}(\Gamma_{D})}\right)^{2}\right]^{\frac{1}{2}}$ (113)
$\displaystyle\leq\left(\lvert
u\rvert^{2}_{H^{1}(\Omega)}+h^{-\frac{\alpha}{2}}C_{a}^{\frac{\alpha}{2}}\lVert\partial_{n}u\rVert_{L^{2}(\Gamma_{D})}+2h^{\frac{\alpha}{2}}C_{a}^{-\frac{\alpha}{2}}\lVert
u\rVert_{L^{2}(\Gamma_{D})}\right)^{\frac{1}{2}}$ (114)
$\displaystyle\cdot\;\;\left(\lvert
v\rvert^{2}_{H^{1}(\Omega)}+h^{-\frac{\alpha}{2}}C_{a}^{\frac{\alpha}{2}}\lVert\partial_{n}v\rVert_{L^{2}(\Gamma_{D})}+2h^{\frac{\alpha}{2}}C_{a}^{-\frac{\alpha}{2}}\lVert
v\rVert_{L^{2}(\Gamma_{D})}\right)^{\frac{1}{2}}$ (115) $\displaystyle\leq
C_{4}\lVert u\rVert_{h,\alpha}\lVert v\rVert_{h,\alpha},$ (116)

hence $C_{4}=2$ is the continuity constant. ∎

######  Corollary 6.2.1.

Using Lax–Milgram lemma we know that for any $(f,g)\in L^{2}(\Omega)\times
H^{\frac{1}{2}}(\partial\Omega)$ there exists a unique $u_{h}\in V_{h}$
solving (83).

######  Proposition 7.

The discrete solution $u_{h}\in V_{h}$ of (83) is a $\delta$-order
approximation of the solution $u\in H^{1+\delta}(\Omega)$ of (82), with
respect to $\lVert\cdot\rVert_{h,\alpha}$ i.e.

$\lVert u-u_{h}\rVert_{h,\alpha}\leq\underset{v_{h}\in V_{h}}{\inf}\lVert
u-v_{h}\rVert_{h,\alpha}\leq C_{5}h^{\delta}\lvert
u\rvert_{H^{1+\delta}(\Omega)},$ (117)

provided (88) holds, $\delta\in(\frac{1}{2},1]$ and
$\alpha\leq\gamma-\beta-1$.

###### Proof.

Via triangle inequality we know that for any $v_{h}\in V_{h}$,

$\displaystyle\lVert u-u_{h}\rVert_{h,\alpha}\leq\lVert
u-v_{h}\rVert_{h,\alpha}+\lVert v_{h}-u_{h}\rVert_{h,\alpha},$ (118)

we can now bound the second term in the inequality, i.e.

$\displaystyle\lVert v_{h}-u_{h}\rVert_{h,\alpha}^{2}\leq
5a_{h}(v_{h}-u_{h},v_{h}-u_{h})=5a_{h}(u-v_{h},u_{h}-v_{h}),$ (119)

where the last identity is a consequence of Galerkin orthogonality, i.e.

$\displaystyle a_{h}(u_{h}-v_{h},u_{h}-v_{h})$
$\displaystyle=a_{h}(u_{h},u_{h}-v_{h})-a_{h}(v_{h},u_{h}-v_{h})$ (120)
$\displaystyle=(f,u_{h}-v_{h})_{L^{2}(\Omega)}-a_{h}(v_{h},u_{h}-v_{h})$ (121)
$\displaystyle=a_{h}(u,u_{h}-v_{h})-a_{h}(v_{h},u_{h}-v_{h})$ (122)
$\displaystyle=a_{h}(u-v_{h},u_{h}-v_{h}).$ (123)

Since $1+\delta>3/2$, from (119) and the continuity of $a_{h}$, we obtain

$\displaystyle\lVert v_{h}-u_{h}\rVert_{h,\alpha}^{2}\leq
5a_{h}(u-v_{h},u_{h}-v_{h})=10\lVert v_{h}-u_{h}\rVert_{h,\alpha}\lVert
u-v_{h}\rVert_{h,\alpha}.$ (124)

Hence we have that $\forall v_{h}\in V_{h}$ the following inequality holds,

$\displaystyle\lVert u-u_{h}\rVert_{h,\alpha}\leq 11\lVert
u-v_{h}\rVert_{h,\alpha}$ (125)

and, since the previous inequality holds for all $v_{h}$ in $V_{h}$ we have
that,

$\displaystyle\lVert u-u_{h}\rVert_{h,\alpha}\leq 11\underset{v_{h}\in
V_{h}}{\inf}\lVert u-v_{h}\rVert_{h,\alpha}.$ (126)

The estimate of the best approximation error for $u\in H^{1+\delta}(\Omega)$
by discrete functions in $V_{h}$, with respect to the
$\lVert\cdot\rVert_{h,\alpha}$ follows, for instance, from the following bound
proved in [31]

$\displaystyle\lVert u-I_{h}u\rVert_{h,\alpha}\leq C_{5}h^{\delta}\lvert
u\rvert_{H^{2}(\Omega)},$ (127)

where $I_{h}:H^{1+\delta}(\Omega)\to V_{h}$ is a special type of Clément
interpolant introduced in [31] and we assume (88) holds. ∎

We can now proceed to use the Aubin–Nitsche duality argument to give an
estimate of the convergence rate of the scheme here proposed with respect to
the $L^{2}(\Omega)$ norm.

######  Proposition 8.

Under the assumption that (88) holds the discrete solution $u_{h}\in V_{h}$ of
(83) is a $2\delta$-order approximation of the solution $u\in
H^{1+\delta}(\Omega)$ of (82) with respect to the $L^{2}(\Omega)$ norm i.e.

$\lVert u-u_{h}\rVert_{L^{2}(\Omega)}\leq C_{6}h^{2\delta}\lvert
u\rvert_{H^{1+\delta}(\Omega)},$ (128)

provided that $\delta\in(\frac{1}{2},1]$ and $\alpha\leq\gamma-\beta-1$.

###### Proof.

We consider the auxiliary problem, find $w\in H^{1}_{\Gamma_{D}}(\Omega)$ such
that

$a(w,z)=(u-u_{h},z)_{L^{2}(\Omega)},\qquad\forall z\in
H^{1}_{\Gamma_{D}}(\Omega).$ (129)

We now notice that

$\displaystyle\lVert u-u_{h}\rVert^{2}_{L^{2}(\Omega)}$
$\displaystyle=\left(\nabla
w,\nabla(u-u_{h})\right)_{L^{2}(\Omega)}-(u-u_{h},\partial_{n}w)_{L^{2}(\partial\Omega)}$
(130) $\displaystyle=\left(\nabla
w,\nabla(u-u_{h})\right)_{L^{2}(\Omega)}-(u-u_{h},\partial_{n}w)_{L^{2}(\Gamma_{D})}$
(131)
$\displaystyle-(\partial_{n}u-\partial_{n}u_{h},w)_{L^{2}(\Gamma_{D})}+C_{a}h^{-1}\left(u-u_{h},w\right)_{L^{2}(\Gamma_{D})}$
(132) $\displaystyle=a_{h}(u-u_{h},w).$ (133)

We now notice that
$a_{h}(u-u_{h},w_{h})=(f,w_{h})_{L^{2}(\Omega)}-(f,w_{h})_{L^{2}(\Omega)}=0$
therefore using the continuity of $a_{h}$ we can rewrite the previous identity
as

$\displaystyle\lVert u-u_{h}\rVert^{2}_{L^{2}(\Omega)}$
$\displaystyle=a_{h}(u-u_{h},w-w_{h})\leq\lVert u-u_{h}\rVert_{h,1}\lVert
w-w_{h}\rVert_{h,1}.$ (134)

Using Proposition 7 we can bound the $L^{2}(\Omega)$ error as

$\displaystyle\lVert u-u_{h}\rVert^{2}_{L^{2}(\Omega)}$ $\displaystyle\leq
C_{5}^{2}h^{\delta}\lvert u\rvert_{H^{2}(\Omega)}h^{\delta}\lvert
w\rvert_{H^{2}(\Omega)}\leq C_{5}h^{2\delta}\lVert
u-u_{h}\rVert_{L^{2}(\Omega)}\lvert u\rvert_{H^{2}(\Omega)}.$ (135)

Dividing the previous inequality by $\lVert u-u_{h}\rVert_{L^{2}(\Omega)}$ we
obtain the desired result. ∎

We now focus our attention to the analysis of the snapping back to the grid
perturbation discussed in Section 5

When adopting a snapping back to grid technique we are solving (83) on a
domain $\widetilde{\Omega}$, which is a perturbation of the original $\Omega$.
We introduce this perturbation to get rid of cells that do not respect the
geometrical assumption (88).

We aim to use the theory presented in [41], and to develop an a priori error
estimate for the snapping back to grid approach viewed as a perturbation
method.

Before diving in the results presented in [41] we need to discuss eventual
regularity issue that might arise from considering the perturbed domain
$\widetilde{\Omega}$. In the one-dimensional case we know that the snapping
presented in Section 5 does not affect the elliptic regularity discussed at
the beginning of this section. Instead, in the two-dimensional case the worst
possible scenario comes when, as a consequence of the snapping back to grid,
the perturbed domain $\widetilde{\Omega}$ presents a re-entrant corner of
aperture $\frac{3}{2}\pi$. An example of such situation can be seen in Figure
4 (a).

In this case elliptic regularity only guarantees that the $\widetilde{u}\in
H^{1}_{\Gamma_{D}}(\Omega)$ solving

$(\nabla\widetilde{u},\nabla
v)_{L^{2}(\widetilde{\Omega})}=(f,v)_{L^{2}(\Omega)}+\prescript{}{H^{-\frac{1}{2}}(\Gamma_{N})}{\langle}g_{N},v\rangle_{H^{\frac{1}{2}}(\Gamma_{N})},\qquad\forall
v\in H^{1}_{\Gamma_{D}}(\Omega),$ (136)

lives in the Sobolev space $H^{\frac{5}{3}+\varepsilon}(\Omega)$ when
$(f,g_{D},g_{N})\in L^{2}(\Omega)\times H^{\frac{1}{2}}(\partial\Omega)\times
H^{-\frac{1}{2}}(\partial\Omega)$ [30]. Hence we might have a mismatch between
the regularity of $\widetilde{u}$ and $u$. For this reason we will consider
the regularity of $\widetilde{u}$ when using estimates of the form of (128).

We can now use the following simplified result from [41, Theorem 8.4] to bound
the $L^{2}$ distance over $R$ between $\widetilde{u}$ and $u$.

######  Theorem 8.1.

Let $\Omega\subset R$ and $\widetilde{\Omega}\subset R$ be Lipschitz domain
with corners of aperture at most $\frac{3}{2}\pi$. We consider $u\in
H^{1}(\Omega)$ and $\widetilde{u}\in H^{1}(\widetilde{\Omega})$ solutions of
(82) and (136) for a data $f\in L^{2}(R)$ and with homogeneous Dirichlet
boundary condition, i.e. $g_{D}=0$ on the entire boundary, i.e.
$\Gamma_{D}=\partial\Omega$ and
$\widetilde{\Gamma}_{D}=\partial\widetilde{\Omega}$ respectively. Under these
hypotheses the following bound for the $L^{2}$ distance over $R$ between
$\widetilde{u}$ and $u$, holds:

$\lVert u-\widetilde{u}\rVert_{L^{2}(R)}\leq
C_{7}\left(\sup_{\mathbf{x}\in\Omega\Delta\widetilde{\Omega}}d\left(\mathbf{x},\partial\Omega\right)\right)\lVert
f\rVert_{L^{2}(R)},$ (137)

where $\Omega\Delta\widetilde{\Omega}$ denotes the set difference between
$\Omega$ and $\widetilde{\Omega}$ and $C_{7}$ is a constant depending only the
diameter of $R$ and the aperture of all the re-entrant corners of
$\widetilde{\Omega}$.

######  Corollary 8.1.1.

Let $\widetilde{u}_{h}\in V_{h}$ be the solution of the perturbed discrete
variational problem,

$\displaystyle(\nabla\widetilde{u}_{h},\nabla{v}_{h})_{L^{2}(\widetilde{\Omega})}+C_{a}h^{-\alpha}(\widetilde{u}_{h},{v}_{h})_{L^{2}(\Gamma_{D})}$
$\displaystyle-(\partial_{n}\widetilde{u}_{h},{v}_{h})_{L^{2}(\Gamma_{D})}-(\widetilde{u}_{h},\partial_{n}{v}_{h})_{L^{2}(\Gamma_{D})}$
(138)
$\displaystyle=(f,{v}_{h})_{L^{2}(\widetilde{\Omega})},\qquad\forall{v}_{h}\in
V_{h},$

where $V_{h}$ is the usual space of bilinear function over each element of a
mesh $\mathcal{T}_{h}$ for $R\subset\mathbb{R}$. If $\widetilde{\Omega}$ is
the geometry obtained from $\Omega$ after snapping back to cell for every $K$
such that

$\lvert\partial\Omega\cap K\rvert\leq\widehat{C}_{\Omega}h^{\beta},$ (139)

then there exists $\delta\in(\frac{2}{3},1+\beta-\gamma]$, depending only on
the regularity of $\partial\widetilde{\Omega}$, such that

$\lVert u-\widetilde{u}_{h}\rVert_{L^{2}(R)}\leq C_{8}h^{2\delta}\lVert
f\rVert_{L^{2}(R)}.$ (140)

###### Proof.

The proof of this result is almost immediate, in fact since we are snapping
back to the grid only when

$\lvert\partial\Omega\cap K\rvert\leq\widehat{C}_{\Omega}h^{\beta},$ (141)

we also know that the following bound holds:

$\sup_{\mathbf{x}\in\Omega\Delta\widetilde{\Omega}}d\left(\mathbf{x},\partial\Omega\right)\leq\sqrt{2}\widehat{C}_{\Omega}h^{\beta}.$
(142)

We then observe that as discussed previously $\widetilde{\Omega}$ can only
have re-entrant corners of aperture $\frac{3}{2}\pi$ and therefore we can
apply the previous theorem to obtain,

$\lVert
u-\widetilde{u}\rVert_{L^{2}(R)}\leq\sqrt{2}\widehat{C}_{\Omega}C_{7}h^{\beta}\lVert
f\rVert_{L^{2}(R)}.$ (143)

Since we know that $\widetilde{u}\in H^{1+\delta}$, with
$\delta\in(\frac{2}{3},1]$, and that (88) is verified by $\mathcal{T}_{h}$ and
$\widetilde{\Omega}$ if we pick $\alpha\leq\gamma-\beta-1$ we can apply (128)
to obtain the following bound,

$\lVert\widetilde{u}-\widetilde{u}_{h}\rVert_{L^{2}(\widetilde{\Omega})}\leq
C_{6}h^{2\delta}\lVert f\rVert_{L^{2}(\widetilde{\Omega})}.$ (144)

Using the fact that
$\lVert\cdot\rVert_{L^{2}(\widetilde{\Omega})}\leq\lVert\cdot\rVert_{L^{2}(R)}$
and the triangular inequality we get that,

$\lVert u-\widetilde{u}_{h}\rVert_{L^{2}(R)}\leq C_{8}h^{2\delta}\lVert
f\rVert_{L^{2}(R)},$ (145)

where $u$ and $\widehat{u}$ are extended to zero outside of $\Omega$ and
$\widehat{\Omega}$ and $C_{8}=(C_{6}+\sqrt{2}\widehat{C}_{\Omega}C_{7})$. ∎

##  9 Numerical tests

In this section we provide some numerical tests showing second order accuracy
in 1D and 2D. We test different smooth geometries: circle, leaf and flower-
shaped domain. Then, we switch to non-smooth geometries, such as hourglass-
shaped domain. For each geometry we show the accuracy of the solution and of
its gradient, that are computed as follows:

$\displaystyle{\rm error}$ $\displaystyle=\frac{||u_{h}-u_{\rm
exa}||_{L^{2}(\mathcal{A})}}{||u_{\rm exa}||_{L^{2}(\mathcal{A})}}$ (146)
$\displaystyle\nabla{\rm error}$ $\displaystyle=\frac{||\nabla u_{h}-\nabla
u_{\rm exa}||_{L^{2}(\mathcal{A})}}{||\nabla u_{\rm
exa}||_{L^{2}(\mathcal{A})}}.$ (147)

In the accuracy analysis, the penalisation parameter plays an important role.
We find a compromise between a second order accurate solution and a well
conditioned numerical scheme.

###  9.1 Numerical tests in 1D

In this section we perform numerical tests to check that the method achieves
the expected accuracy independently of the parameters $\theta_{1}$ and
$\theta_{2}$ (see Fi. 1).

We choose the following exact solution

$u_{\rm exa}=\sin(k_{0}x+k_{1}),\qquad{\rm in}\,[0,1]$ (148)

with $k_{0}=5$, $k_{1}=1$, so that both Dirichlet and Neumann conditions are
non zero at the boundaries. In Fig. 6, we show the comparison between the
exact solution (blue line) and the numerical solution (red asterisks) for the
mixed problem in Eq. (38), with Dirichlet boundary conditions in $a$, Neumann
boundary conditions in $b$, and number of elements $N=20$. In panels (b), (c),
(d) we show the convergence order. while in (e), (f), (g) the slope of the
conditioning number of the stiff matrix for a fixed value of
$\theta_{2}=10^{-3}$,
$\theta_{1}\in\\{0.0010,0.0015,0.0020,0.0025,\cdots,0.990\\}$ and
$N\in\\{20,40,80,160,320,640\\}$. In this test the penalisation coefficient is
$\lambda=h^{-\alpha}$, where, in panels (b),(e) $\alpha=2$, in (c),(f)
$\alpha=1.75$ and in (d),(g) $\alpha=1.5$.

\begin{overpic}[abs,width=198.7425pt,unit=1mm,scale={.25}]{Figures/solu_sin_1D.png}
\put(-2.0,48.0){(a)}\end{overpic}

\begin{overpic}[abs,width=437.2359pt,unit=1mm,scale={.25}]{Figures/accuracy_order_fixed_theta2_alpha2}
\put(0.0,34.0){(b)} \end{overpic}

\begin{overpic}[abs,width=437.2359pt,unit=1mm,scale={.25}]{Figures/accuracy_order_fixed_theta2_alpha1_75}
\put(0.0,34.0){(c)} \end{overpic}

\begin{overpic}[abs,width=437.2359pt,unit=1mm,scale={.25}]{Figures/accuracy_order_fixed_theta2_alpha1_5}
\put(0.0,34.0){(d)} \end{overpic}

\begin{overpic}[abs,width=437.2359pt,unit=1mm,scale={.25}]{Figures/condA_order_fixed_theta2_alpha2}
\put(0.0,34.0){(e)} \end{overpic}

\begin{overpic}[abs,width=437.2359pt,unit=1mm,scale={.25}]{Figures/condA_order_fixed_theta2_alpha1_75}
\put(0.0,34.0){(f)} \end{overpic}

\begin{overpic}[abs,width=437.2359pt,unit=1mm,scale={.25}]{Figures/condA_order_fixed_theta2_alpha1_5}
\put(0.0,34.0){(g)} \end{overpic}

Figure 6: (a): Numerical solution of the mixed problem (62) (asterisks) and
comparison with the exact solution (148), with $k_{0}=5.0$ and $k_{1}=1.0$.
$N=20$ intervals have been adopted. Convergence order (b),(c),(d) and
conditioning slope (e),(f), (g) for a fixed value of $\theta_{2}=10^{-3}$,
$\theta_{1}\in\\{0.0010,0.0015,0.0020,0.0025,\cdots,0.990\\}$ and
$N\in\\{20,40,80,160,320,640\\}$. In this test $\lambda=h^{-\alpha}$, where,
in panels (b),(e) $\alpha=2$, in (c),(f) $\alpha=1.75$ and in (d),(g)
$\alpha=1.5$.

###  9.2 Numerical results in 2D

In this section we provide some numerical tests showing second order accuracy
in two dimensions. We test different smooth geometries: circle, flower and
leaf-shaped domain, and, at the end, a saddle point domain. The rectangular
region $R\subset\mathbb{R}^{2}$, and the different domains $\Omega\subset R$,
are defined by level-set functions that we define below. In Figs. 7-11, we
show, in panel (a), the shape of the domain that we are considering, in panel
(b) the accuracy order of the solution for different values of penalisation
coefficient $\lambda=h^{-\alpha}$ and snapping parameters $1/\lambda$ (defined
in Algorithm 2), with $\alpha\in\\{1.5,1.75,2\\}$. Furthermore, in panel (c),
we show the slope of the conditioning number of the stiff matrix, and in panel
(d), the accuracy order of the gradient of the solution.

#### Circular domain

Let us start with a circular domain $\Omega\subset[0,1]^{2}$ centered at
$(x_{c},y_{c})\in\Omega$. To prove that the order of accuracy does not depend
on the symmetry of the domain, we calculate 10 different numerical solutions
coming from different level-set functions. For this reason, we choose the
following expression for the center of the domain:

$(x_{c},y_{c})=\left(0.5+\widehat{\epsilon}_{i}\Delta
x,0.5+\widetilde{\epsilon}_{j}\Delta y\right),\quad i\in 1,\cdots,10,$

where $\widehat{\epsilon}_{i},\widetilde{\epsilon}_{j}\in[0,1]$, are randomly
chosen. In this case, the level-set function is
$\phi=r-\sqrt{(x-x_{c})^{2}+(y-y_{c})^{2}}$, where $r=0.4$ is the radius of
the circle. At the end we have 10 different variables for the quantities
defined in Eqs. (146-147), ${\rm error},\nabla{\rm error}$, one for each pair
of $(x_{c},y_{c})$. In Figs. 7-8 we show the average of these 10 quantities.

The exact solution we choose is $u_{\rm exa}=\cos(2\pi x)\cos(2\pi y)$ with
homogeneous Dirichlet boundary conditions (i.e., $\Gamma_{D}=\Gamma$) in Fig.
7, and mixed boundary conditions in Fig. 8, i.e.
$\Gamma_{D}=\Gamma\cap\\{x\leq 0.5\\}$ and $\Gamma_{N}=\Gamma\cap\\{x>0.5\\}$.

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/circle_Domain.png}
\put(-3.0,48.0){(a)} \end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_snapping_halpha.png}
\put(-2.0,48.0){(b)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/condA_snapping_halpha.png}
\put(-2.0,48.0){(c)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_gradient_snapping_halpha.png}
\put(-2.0,48.0){(d)}\end{overpic}

Figure 7: Numerical results for the circular domain and Dirichlet boundary
conditions: (a) shape of the domain; (b) error plot for different values of
the snapping back to grid and penalisation term; (c) conditioning number for
different values of the snapping back to grid and penalisation term; (d) error
plot for the gradient of the solution for different values of the snapping
back to grid and penalisation term.

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/difference_solu_circle_mixed.png}
\put(-3.0,48.0){(a)} \end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_snapping_alpha_circle_mixed.png}
\put(-2.0,48.0){(b)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/condA_snapping_halpha_circle_mixed.png}
\put(-2.0,48.0){(c)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_grad_alpha_circle_mixed.png}
\put(-2.0,48.0){(d)}\end{overpic}

Figure 8: Numerical results for the circular domain and mixed boundary
conditions: (a) difference $|u_{h}-u_{\rm exa}|$; (b) error plot for different
values of the snapping back to grid and penalisation term; (c) conditioning
number for different values of the snapping back to grid and penalisation
term; (d) error plot for the gradient of the solution for different values of
the snapping back to grid and penalisation term.

#### Flower-shaped domain

Here we show a flower-shaped domain $\Omega\subset R=[-1,1]^{2}$, the exact
solution we choose is $u_{\rm exa}=\exp(xy)$ and homogeneous Dirichlet
boundary conditions, i.e. $\Gamma_{D}=\Gamma$. In this case, the level-set
function is

$\displaystyle X=x-0.03\sqrt{3},\quad Y=y-0.04\sqrt{2},\quad
R=\sqrt{X^{2}+Y^{2}}$
$\displaystyle\phi=\frac{R-0.52-(Y^{5}+5X^{4}Y-10X^{2}Y^{3})}{5R^{5}},$

and the results are showed in Fig. 9.

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/flower_Domain.png}
\put(-3.0,48.0){(a)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_snapping_halpha_flower.png}
\put(-2.0,48.0){(b)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/condA_snapping_halpha_flower.png}
\put(-2.0,48.0){(c)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_gradient_flower_snapping_halpha.png}
\put(-2.0,48.0){(d)}\end{overpic}

Figure 9: Numerical results for the flower domain: (a) shape of the domain;
(b) error plot for different values of the snapping back to grid and
penalisation term; (c) conditioning number for different values of the
snapping back to grid and penalisation term; (d) error plot for the gradient
of the solution for different values of the snapping back to grid and
penalisation term.

#### Leaf domain

Here we show a leaf domain $\Omega\subset R=[0,1]^{2}$, the exact solution we
choose is $u_{\rm exa}=\exp(xy)$ and mixed boundary conditions, such that
$\Gamma_{D}=\partial\Omega\cap\\{x<0.5\\}$ and
$\Gamma_{N}=\partial\Omega\cap\\{x\geq 0.5\\}$. In this case, the level-set
function is

$\displaystyle x_{1}=0.4,\quad y_{1}=0.5,\quad
R_{1}=\sqrt{(y-x_{1})^{2}+(y-y_{1})^{2}}$ $\displaystyle x_{2}=0.6,\quad
y_{2}=0.5,\quad R_{2}=\sqrt{(x-x_{2})^{2}+(y-y_{2})^{2}}$ $\displaystyle
r_{0}=0.4,\quad\phi_{1}=R_{1}-r_{0},\quad\phi_{2}=R_{2}-r_{0},\quad\phi=\max\\{\phi_{1},\phi_{2}\\},$

and the results are in Fig. 10.

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/leaf_Domain.png}
\put(-3.0,48.0){(a)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_snapping_alpha_leaf.png}
\put(-2.0,48.0){(b)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/condA_snapping_halpha_leaf.png}
\put(-2.0,48.0){(c)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_gradient_leaf_snapping_halpha.png}
\put(-2.0,48.0){(d)}\end{overpic}

Figure 10: Numerical results for the leaf domain: (a) shape of the domain; (b)
error plot for different values of the snapping back to grid and penalisation
term; (c) conditioning number for different values of the snapping back to
grid and penalisation term; (d) error plot for the gradient of the solution
for different values of the snapping back to grid and penalisation term.

#### Hourglass domain

Here we consider a domain $\Omega\subset R=[-1,1]^{2}$, with a saddle point,
whose level-set function is

$\displaystyle X=x-0.03\sqrt{3},\qquad
Y=y-0.04\sqrt{2},\qquad\phi(x,y)=256\,Y^{4}-16\,X^{4}-128\,Y^{2}+36\,X^{2}.$

The exact solution we choose is $u_{\rm exa}=\exp(xy)$ with homogeneous
Dirichlet boundary conditions (i.e., $\Gamma_{D}=\Gamma$) in Fig. 11, and
mixed boundary conditions in Fig. 12, i.e. $\Gamma_{D}=\Gamma\cap\\{x\leq
0\\}$ and $\Gamma_{N}=\Gamma\cap\\{x>0\\}$.

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/hourglass_Domain.png}
\put(-3.0,48.0){(a)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_snapping_alpha_hourglass.png}
\put(-2.0,48.0){(b)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/condA_snapping_halpha_hourglass.png}
\put(-2.0,48.0){(c)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_gradient_hourglass_snapping_halpha.png}
\put(-2.0,48.0){(d)}\end{overpic}

Figure 11: Numerical results for the hourglass domain: (a) shape of the
domain; (b) error plot for different values of the snapping back to grid and
penalisation term; (c) conditioning number for different values of the
snapping back to grid and penalisation term; (d) error plot for the gradient
of the solution for different values of the snapping back to grid and
penalisation term.

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/difference_solu_hourglass_mixed.png}
\put(-3.0,48.0){(a)} \end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_snapping_alpha_hourglass_mixed.png}
\put(-2.0,48.0){(b)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/cA_snapping_alpha_hourglass_mixed.png}
\put(-2.0,48.0){(c)}\end{overpic}

\begin{overpic}[abs,width=397.48499pt,unit=1mm,scale={.25}]{Figures/error_grad_snapping_alpha_hourglass_mixed.png}
\put(-2.0,48.0){(d)}\end{overpic}

Figure 12: Numerical results for the hourglass domain and mixed boundary
conditions: (a) difference $|u_{h}-u_{\rm exa}|$; (b) error plot for different
values of the snapping back to grid and penalisation term; (c) conditioning
number for different values of the snapping back to grid and penalisation
term; (d) error plot for the gradient of the solution for different values of
the snapping back to grid and penalisation term.

##  10 Conclusions

In the present work we have presented a novel nodal ghost discretisation posed
on regular square meshes. In particular, our derivation of the numerical
scheme here presented found its root in the variational formulation of the
Poisson equation. Such formulation, often known also as weak formulation, is
at heart of the finite element method.

Resorting to a variational formulation allows us to preserve the symmetry and
the coercivity of the continuous problem in the discrete setting. The
conservation of such properties at discrete level represents a substantial
difference and a great improvement with respect to other nodal ghost schemes
such as the Coco-Russo method [19, 20].

There are many difference between the method here proposed and a standard
finite element discretisation. One particular difference, between our scheme
and a standard body fitted FEM, is that we described the geometry of the
domain $\Omega$ via a level-set function. Such approach is also adopted in
CutFEM, [31, 14, 15, 16], which differs from our scheme for the following
reasons. First of all, in the CutFEM in order to compute the boundary
integrals appearing, they introduce a sub-triangulation of cut elements. Such
process requires additional computational work. We avoid the construction of a
sub-triangulation via the quadrature technique presented in Section 5.
Secondly, we escape the so called small cut issue, which might result in
losing the coercivity of the bilinear form, by perturbing mesh geometry. In
particular, the perturbation procedure we adopted is well known to engineering
community by the name of snapping back to grid.

Adopting standard techniques, in the FEM and CutFEM, we are able to prove
optimal converge rate for our scheme with respect to the $L^{2}$ norm. To the
best of the authors knowledge up until now, there has not been a rigorous
study of a priori convergence rates when the snapping back to grid technique
is adopted. For this reason, we have presented a novel approach to prove a
priori convergence rates for our scheme even when the snapping back to grid
technique is used to address the small cuts issue.

Lastly, we have presented many numerical experiments that validate the $L^{2}$
convergence rate theoretically proved. Furthermore, we would like to point
that the numerical experiments shown in Section 9 suggest that the error in
energy norm is super-converging with respect to the above mentioned analysis.

In future, we plan to further investigate the scheme here proposed. In
particular, we would like to study in greater details the eigenvalue problem
(84) and the role it plays when determining the optimal value of the
penalisation constant. We are also interested in developing fast solvers for
our discretisation based on a multigrid approach, similarly to [19].

## References

  * [1] Clarissa Astuto, Armando Coco, and Giovanni Russo.  A finite-difference ghost-point multigrid method for multi-scale modelling of sorption kinetics of a surfactant past an oscillating bubble.  _Journal of Computational Physics_ , 476:111880, 2023.
  * [2] Clarissa Astuto, Mohammed Lemou, and Giovanni Russo.  Time multiscale modeling of sorption kinetics i: uniformly accurate schemes for highly oscillatory advection-diffusion equation.  _arXiv preprint arXiv:2307.14001_ , 2023.
  * [3] Clarissa Astuto,  High order multiscale methods for advection-diffusion equation with highly oscillatory boundary condition.  _arXiv preprint arXiv:2401.12226_ , 2024\.
  * [4] Clarissa Astuto, Antonio Raudino, and Giovanni Russo.  Multiscale modeling of sorption kinetics.  _Multiscale Modeling & Simulation_, 21(1):374–399, 2023.
  * [5] Jean-Pierre Aubin.  Approximation of elliptic boundary-value problems.  _Wiley-Interscience [A division of John Wiley & Sons, Inc.]_, New York-London-Sydney, 1972.  _Pure and Applied Mathematics, Vol. XXVI_.
  * [6] Ivo Babuška.  The finite element method with penalty.  _Math. Comp._ , 27:221–228, 1973.
  * [7] John W. Barrett and Charles M. Elliott.  Finite element approximation of the Dirichlet problem using the boundary penalty method.  _Numer. Math._ , 49(4):343–366, 1986.
  * [8] Daniele Boffi, Andrea Cangiani, Marco Feder, Lucia Gastaldi, and Luca Heltai.  A comparison of non-matching techniques for the finite element approximation of interface problems.  _Computers & Mathematics with Applications_, 151:101–115, 2023.
  * [9] Daniele Boffi and Lucia Gastaldi.  A fictitious domain approach with Lagrange multiplier for fluid-structure interactions.  _Numerische Mathematik_ , 135(3):711–732, 2017.
  * [10] Daniele Boffi, Lucia Gastaldi, Luca Heltai, and Charles S. Peskin.  On the hyper-elastic formulation of the immersed boundary method.  _Computer Methods in Applied Mechanics and Engineering_ , 197(25):2210–2231, 2008.  _Immersed Boundary Method and Its Extensions_.
  * [11] Susanne C. Brenner and L. Ridgway Scott.  _The Mathematical Theory of Finite Element Methods_ , volume 15 of Texts in Applied Mathematics.  _Springer New York_ , New York, NY, 2008.
  * [12] Haim Brezis.  _Functional Analysis, Sobolev Spaces and Partial Differential Equations_ , Universitext.  _Springer Nature_ , New York, NY, 1 edition, 2010.
  * [13] Th Bulenda and Jan Knippers.  Stability of grid shells.  _Computers & Structures_, 79(12):1161–1174, 2001.
  * [14] Erik Burman, Susanne Claus, Peter Hansbo, Mats G. Larson, and André Massing.  Cutfem: Discretizing geometry and partial differential equations.  _International Journal for Numerical Methods in Engineering_ , 104(7):472–501, 2015.
  * [15] Erik Burman and Peter Hansbo.  Fictitious domain finite element methods using cut elements: II. A stabilized Nitsche method.  Applied Numerical Mathematics, 62(4):328-341, 2012. Third Chilean Workshop on Numerical Analysis of Partial Differential Equations (WONAPDE 2010).
  * [16] Erik Burman, Peter Hansbo, and Mats G. Larson.  Cutfem based on extended finite element spaces.  _Numerische Mathematik_ , 152(2):331–369, Oct 2022.
  * [17] Han Chen, Chohong Min, and Frédéric Gibou.  A supra-convergent finite difference scheme for the Poisson and heat equations on irregular domains and non-graded adaptive cartesian grids.  _Journal of scientific computing_ , 31(1-2):19–60, 2007.
  * [18] Philippe G. Ciarlet.  _The finite element method for elliptic problems_.  _Classics in applied mathematics ; 40. Society for Industrial and Applied Mathematics SIAM_ , 3600 Market Street, Floor 6, Philadelphia, PA 19104, Philadelphia, Pa, 2002.
  * [19] Armando Coco and Giovanni Russo.  Finite-difference ghost-point multigrid methods on cartesian grids for elliptic problems in arbitrary domains.  _Journal of computational physics_ , 241:464–501, 2013.
  * [20] Armando Coco and Giovanni Russo.  Second order finite-difference ghost-point multigrid methods for elliptic problems with discontinuous coefficients on an arbitrary interface.  _Journal of computational physics_ , 361:299–330, 2018.
  * [21] Armando Coco.  A multigrid ghost-point level-set method for incompressible Navier-Stokes equations on moving domains with curved boundaries.  _Journal of computational physics_ , 2020.
  * [22] R. Courant.  Variational methods for the solution of problems of equilibrium and vibrations.  _Bull. Amer. Math. Soc._ , 49:1–23, 1943.
  * [23] R. Courant and D. Hilbert.  _Methods of mathematical physics. Vol. I_.  _Interscience Publishers, Inc., New York, N.Y._ , 1953.
  * [24] L. Beirão da Veiga, A. Buffa, G. Sangalli, and R. Vázquez.  Mathematical analysis of variational isogeometric methods.  _Acta Numerica_ , 23:157–287, 2014.
  * [25] Herbert Edelsbrunner.  Triangulations and meshes in computational geometry.  _Acta Numerica_ , 9:133–213, 2000.
  * [26] Herbert Edelsbrunner.  _Geometry and Topology for Mesh Generation_.  _Cambridge Monographs on Applied and Computational Mathematics. Cambridge University Press_ , 2001.
  * [27] Alexandre Ern and Jean-Luc Guermond.  _Finite Elements I: Approximation and Interpolation_ , volume 72 of Texts in Applied Mathematics.  _Springer International Publishing, Cham_ , 2021.
  * [28] Frederic Gibou, Ronald P. Fedkiw, Li-Tien Cheng, and Myungjoo Kang.  A second-order-accurate symmetric discretization of the Poisson equation on irregular domains.  _Journal of computational physics_ , 176(1):205–227, 2002.
  * [29] Roland Glowinski, Tsorng-Whay Pan, and Jacques Periaux.  A fictitious domain method for Dirichlet problem and applications.  _Computer Methods in Applied Mechanics and Engineering_ , 111(3):283–303, 1994.
  * [30] Pierre Grisvard.  _Elliptic problems in nonsmooth domains_.  _SIAM_ , 2011.
  * [31] Anita Hansbo and Peter Hansbo.  An unfitted finite element method, based on Nitsche’s method, for elliptic interface problems.  _Computer Methods in Applied Mechanics and Engineering_ , 191(47):5537–5552, 2002.
  * [32] William D. Henshaw.  Automatic grid generation.  _Acta Numerica_ , 5:121–148, 1996.
  * [33] T.J.R. Hughes, J.A. Cottrell, and Y. Bazilevs.  Isogeometric analysis: CAD, finite elements, NURBS, exact geometry and mesh refinement.  _Computer Methods in Applied Mechanics and Engineering_ , 194(39):4135–4195, 2005.
  * [34] Christoph Lehrenfeld.  High order unfitted finite element methods on level set domains using isoparametric mappings.  _Computer Methods in Applied Mechanics and Engineering_ , 300:716–733, 2016.
  * [35] Christoph Lehrenfeld and Arnold Reusken.  Analysis of a high order unfitted finite element method for an elliptic interface problem.  _IMA J. Numer. Anal._ , 38:1351–1387, 2018.
  * [36] Zi-Cai Li.  Penalty-combined approaches to the Ritz-Galerkin and finite element methods for singularity problems of elliptic equations.  _Numerical Methods for Partial Differential Equations_ , 8(1):33–57, 1992.
  * [37] Amedeo Manuello, Fabio Bazzucchi, and Alberto Carpinteri.  Step-by-step stability analysis of shallow grid shells: Buckling versus snap-through.  In _Proceedings of IASS Annual Symposia_ , volume 2016, pages 1–10. International Association for Shell and Spatial Structures (IASS), 2016.
  * [38] J. 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(1):9–15, July 1971.
  * [39] Stanley Osher and Ronald P. Fedkiw.  _Level set methods and dynamic implicit surfaces_ , Applied mathematical sciences ; 153.  _Springer, New York_ , 2002.
  * [40] Giovanni Russo and Peter Smereka.  A remark on computing distance functions.  _Journal of computational physics_ , 163(1):51–67, 2000.
  * [41] Giuseppe Savaré and Giulio Schimperna.  Domain perturbations and estimates for the solutions of second-order elliptic equations.  _Journal de mathématiques pures et appliquées_ , 81(11):1071–1112, 2002.
  * [42] Joachim Schöberl.  Netgen: an advancing front 2d/3d-mesh generator based on abstract rules.  _Computing and Visualization in Science_ , 1(1):41–52, Jul 1997.
  * [43] J. A. Sethian.  _Level set methods and fast marching methods: evolving interfaces in computational geometry, fluid mechanics, computer vision, and materials science_.  _Cambridge monographs on applied and computational mathematics 3. Cambridge University Press_ , 2nd ed edition, 1999.
  * [44] G. H. Shortley and R. Weller.  The numerical solution of Laplace’s equation.  _Journal of applied physics_ , 9(5):334–348, 1938.
  * [45] Mark Sussman, Peter Smereka, and Stanley Osher.  A level set approach for computing solutions to incompressible two-phase flow.  _Journal of computational physics_ , 114(1):146–159, 1994.