# A variant of the Raviart-Thomas method to handle smooth domains using straight-edged triangles

Fleurianne Bertrand, Vitoriano Ruas

► To cite this version:

Fleurianne Bertrand, Vitoriano Ruas. A variant of the Raviart-Thomas method to handle smooth domains using straight-edged triangles. ESAIM: Mathematical Modelling and Numerical Analysis, 2025, 59, pp.1791-1829. 10.1051/m2an/2025028. hal-05451470

HAL Id: hal-05451470

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

Submitted on 28 Jan 2026

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

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

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

Distributed under a Creative Commons CC BY 4.0 - Attribution - International License
# A VARIANT OF THE RAVIART–THOMAS METHOD TO HANDLE SMOOTH DOMAINS USING STRAIGHT-EDGED TRIANGLES

FLEURIANNE BERTRAND<fnref n="1" /> AND VITORIANO RUAS<fnref n="2" />, *

**Abstract.** Several physical problems modeled by second-order elliptic equations can be efficiently solved using mixed finite elements of the Raviart–Thomas family $RT_k$ for $N$-simplexes, introduced in the seventies. In case Neumann conditions are prescribed on a curvilinear boundary, the normal component of the flux variable should preferably not take up values at nodes shifted to the boundary of the approximating polytope in the corresponding normal direction. This is because the method’s accuracy downgrades, which was shown in Bertrand and Starke (SIAM J. Numer. Anal. 54 (2016) 3648–3667). In that work an order-preserving technique was studied, based on a parametric version of these elements with curved simplexes. In this paper an alternative with straight-edged triangles for two-dimensional problems is proposed. The key point of this method is a Petrov–Galerkin formulation of the mixed problem, in which the test-flux space is a little different from the shape-flux space. After describing the underlying variant of $RT_k$, we show that it gives rise to uniformly stable methods taking the Poisson equation as a model problem. Then a priori error estimates of optimal order in the $[H(div) \times L^2]$-norm are demonstrated.

**Mathematics Subject Classification.** 65N30, 74S05, 76M10, 78M10, 80M10.

Received January 31, 2024. Accepted March 31, 2025.

## 1. INTRODUCTION

This work deals with a particular type of Petrov–Galerkin formulation, primarily designed to preserve the order of finite element methods in the solution of boundary value problems posed in a two- or a three-dimensional domain having a smooth boundary of arbitrary shape, on which DOFs (degrees of freedom) are prescribed. In [27,28] it was applied to conforming Lagrange and Hermite finite element methods with straight-edged triangles and tetrahedra to solve both second and fourth order elliptic equations, in the aim of improving their accuracy, while maintaining orders higher than one in the inherent energy norm. In [31] it was shown that this technique has the very same effect, as applied to nonconforming finite element methods, even those having boundary-prescribed DOFs different from function values at vertexes or any other point of an edge or a face of the elements. Although this approach can be adopted for any type of element geometry, it appears most effective

*Keywords and phrases.* Asymmetric mixed formulations, curved domains, error estimates, finite elements, flux boundary conditions, Poisson equation, Raviart–Thomas, straight-edged triangles.

<footnote n="1">¹ Fakultät für Mathematik, Technische Universität Chemnitz, Chemnitz, Germany.</footnote>

<footnote n="2">² Institut Jean Le Rond d'Alembert, CNRS UMR 7190, Sorbonne Université, Campus Pierre et Marie Curie, Campus Pierre et Marie Curie, 4 place jussieu, Couloir 55-65, 4ème étage, 75005 Paris, France.</footnote>

<footnote>*Corresponding author: vitoriano.ruas@upmc.fr</footnote>
as combined to simplexes. The technique under consideration is simple to implement in both two- and three-
dimensional geometries, and optimal orders can be attained very naturally in both cases. Moreover, in the case
of simplicial elements it definitively eliminates the need for curved elements, and consequently the use of non
affine mappings.

Here we further illustrate the application of this formulation, in the framework of the family of Raviart–Thomas mixed finite elements for simplexes to solve second order boundary value problems in a smooth plane domain, with Neumann conditions prescribed, at least on a boundary portion. We recall that in this case the normal components of the primal (flux) variable in the underlying mixed formulation are prescribed on such a boundary portion. As a matter of fact, with the current proposal we attempt to show once more that the approximate flux variable should preferably not take up corresponding prescribed values at nodes shifted to the boundary of the approximating polytope. Actually, even when there is no node shift, the prescribed flux on the boundary should be normal to the true boundary and not to the approximate boundary. This is because in the former case, a significant increase of accuracy occurs. At this point we should report that a technique based on this principle was studied in [5, 6] for Raviart–Thomas mixed elements. It is based on a parametric version of those elements with curved simplexes. In contrast, in this paper the aforementioned Petrov–Galerkin formulation with straight-edged triangles is employed instead to tackle this problem.

In this article we focus on the uniform stability and the convergence analysis of our method. With this aim we recall that, according to the well known theory set up by Babuška [2, 3] and Brezzi [9], a linear variational problem is well-posed only if the underlying continuous bilinear form is weakly coercive. In the case of problems defined upon a pair of spaces with the same finite dimension, such as finite-element spaces, this reduces to proving that the underlying bilinear form satisfies an *inf-sup* condition for the pair of spaces under consideration (see e.g., [16]). Here, in order to achieve this, we resort to the theoretical framework supplied in [30] for the case of asymmetric mixed formulations. By this denomination we mean mixed formulations in which the shape-product space is different from the test-product space.

In order to bypass non essential difficulties, we study our methodology in connection with the Raviart–Thomas mixed finite elements, taking as a model the Poisson first order system in a bounded plane domain $\Omega$.

An outline of the paper is as follows. In Section 2 we define the model problem to be treated by our mixed finite element method. General considerations are made on the kind of Petrov–Galerkin formulation employed in this work to approximate it and corresponding error bounds are given. This material is followed by the presentation in Section 3 of concepts, notations and properties related to the family of meshes considered in the sequel, more particularly in connection with the boundary of the curved domain. Section 4 is devoted to the description of the specific spaces used here to approximate the model problem with mixed or Neumann boundary conditions, namely, Raviart–Thomas spaces of arbitrary order based on straight-edged triangles; corresponding feasibility results are proven. In Section 5 the well-posedness and uniform stability analysis of the approximate problem is carried out. Section 6 deals with the interpolation errors in the $H(div)$-norm for the space of shape-fluxes under consideration. This leads to best possible a priori error estimates in the natural norm established in Section 7. In Section 8 we supply a relevant complement, make some comments and provide a short summary of the article.

We left aside lengthy proofs for a result used in Section 5, but for the sake of completeness they are given in Appendix I. Additionally we report in Appendix II some numerical experiments confirming formal results obtained in this work.

Throughout our study we shall use the following notations, most related to well known Sobolev function
spaces (see e.g., [1]), extended to the corresponding vector versions: First of all $\mathcal{A} \cdot \mathcal{B}$ represents the inner
product of two tensors $\mathcal{A}$ and $\mathcal{B}$ of order $n$ for $n \in \mathbb{N}$. $D$ being a bounded subset of $\mathbb{R}^n$, $\| \cdot \|_{r,D}$ and $| \cdot |_{r,D}$
denote the standard norm and semi-norm of Sobolev space $H^r(D)$, for $r \in \mathbb{R}^+$ with $H^0(D) = L^2(D)$. We further
denote by $\| \cdot \|_{m,p,D}$ and $| \cdot |_{m,p,D}$ the usual norm and semi-norm of $W^{m,p}(D)$ for $m \in \mathbb{N}^*$ and $p \in [1, \infty] \setminus \{2\}$
with $W^{0,p}(D) = L^p(D)$, and eventually for $W^{m,2}(D) = H^m(D)$ as well. In case $D = \Omega$ we omit this subscript.
$P_k(D)$ represents the space of polynomials of degree less than or equal to $k$ defined in $D$. Finally we introduce
the notation $D^j\mathcal{A}$ for the $j$-th (resp. $(j+1)$-th) order tensor, whose components are the $j$-th order partial derivatives with respect to the space variables of a function (resp. a vector field) $\mathcal{A}$ in the strong or in the weak sense.

## 2. THE MODEL PROBLEM AND ITS ABSTRACT PETROV–GALERKIN APPROXIMATION

In our reliability study we shall apply the abstract results given in [30] to a variant of the Raviart–Thomas family of mixed elements for triangles, specially designed to handle smooth two-dimensional domain. We take the Poisson equation with mixed or Neumann boundary conditions as a model, in the particular framework that we next describe.

Let $\Omega$ be a two-dimensional domain assumed to be at least of the piecewise $C^2$-class and $\Gamma$ be its boundary, consisting of two portions $\Gamma_0$ and $\Gamma_1$ with $length(\Gamma_1) > 0$. We denote by $\mathbf{n}$ the outer normal vector to $\Gamma$ and by $V$, either the space $L^2(\Omega)$ if $length(\Gamma_0) > 0$, or its subspace $L_0^2(\Omega)$ otherwise, where $L_0^2(D) := \{g|g \in L^2(D), \int_D g = 0\}$ for a given domain $D$ of $\mathbb{R}^2$.

We consider the following model equation:

$$ 
\left\{ 
\begin{array}{ll} 
\text{Given } f \in V \text{ find } (\mathbf{p}; u) \text{ such that,} \\ 
-\nabla \cdot \mathbf{p} = f & \text{in } \Omega; \\ 
\mathbf{p} - \nabla u = \mathbf{0} & \text{in } \Omega; \\ 
u = 0 & \text{on } \Gamma_0 \text{ if } length(\Gamma_0) > 0 \text{ or } \int_{\Omega} u = 0 \text{ otherwise;} \\ 
\mathbf{p} \cdot \mathbf{n} = 0 & \text{on } \Gamma_1. 
\end{array} 
\right. 
\quad (1) 
$$

Let $\mathbf{Q}$ be the subspace of $\mathbf{H}(\text{div}, \Omega)$ of those fields $\mathbf{q}$ such that $\mathbf{q} \cdot \mathbf{n} = 0$ on $\Gamma_1$ (cf. [21]). Denoting by $(\cdot, \cdot)$ the standard inner product of $L^2(\Omega)$ in scalar or vector version, problem (1) can be recast in the usual well-posed equivalent variational form (2) and (3) (cf. [25]) that we next recall.

$$ \text{Find } (\mathbf{p}; u) \in \mathbf{Q} \times V \text{ such that } c((\mathbf{p}; u), (\mathbf{q}; v)) = (f, v) \quad \forall (\mathbf{q}; v) \in \mathbf{Q} \times V, \quad (2) $$

where

$$ c((\mathbf{r}; w), (\mathbf{q}; v)) := (\mathbf{r}, \mathbf{q}) + (w, \nabla \cdot \mathbf{q}) - (\nabla \cdot \mathbf{r}, v) \quad \forall ((\mathbf{r}; w); (\mathbf{q}; v)) \in (\mathbf{Q} \times V) \times (\mathbf{Q} \times V). \quad (3) $$

In an attempt to handle the boundary conditions for $\mathbf{p}$ more accurately, in this work we will employ a kind of Petrov–Galerkin finite-element formulation to approximate problem (2) and (3), which we next describe.

To begin with, assume that $\Omega$ is approximated by a polygon $\Omega_h$ defined as the union of triangles in a mesh with maximum edge length equal to $h$. Let us denote by $(\cdot, \cdot)_h$ the standard inner product of $L^2(\Omega_h)$. While $\| \cdot \|_{0,p,h}$ represents the standard norm of $L^p(\Omega_h)$ with $p \in [1, +\infty] \setminus \{2\}$, let $\| \cdot \|_{0,h}$ stand for the norm of $L^2(\Omega_h)$ associated with $(\cdot, \cdot)_h \cdot [\cdot]_h$ being the standard norm of $\mathbf{H}(\text{div}, \Omega_h)$ we denote by $\|\cdot\|_h$ the norm of $\mathbf{H}(\text{div}, \Omega_h) \times L^2(\Omega_h)$ given by $\|\| (\mathbf{q}; v) \|\|_h := [[\mathbf{q}]\|_h + \|v\|_{0,h} \quad \forall (\mathbf{q}; v) \in \mathbf{H}(\text{div}, \Omega_h) \times L^2(\Omega_h)$.

We consider a pair of spaces $(\mathbf{P}_h; \mathbf{Q}_h) \subset [\mathbf{H}(\text{div}; \Omega_h)]^2$ together with a subspace $V_h$ of $L^2(\Omega_h)$ and also of $L_0^2(\Omega_h)$ if $length(\Gamma_0) = 0$, all of them being finite-dimensional.

In our approximation of (2) and (3) the bilinear form $c$ will be replaced by $c_h$ given by

$$ c_h((\mathbf{r}; w), (\mathbf{q}; v)) := a_h(\mathbf{r}, \mathbf{q}) + b_h(w, \mathbf{q}) + d_h(\mathbf{r}, v) \quad \forall (\mathbf{r}; w) \in \mathbf{P}_h \times V_h \text{ and } \forall (\mathbf{q}; v) \in \mathbf{Q}_h \times V_h \quad (4) $$

where

$$ a_h(\mathbf{r}, \mathbf{q}) := (\mathbf{r}, \mathbf{q})_h; \quad b_h(w, \mathbf{q}) := (w, \nabla \cdot \mathbf{q})_h; \quad d_h(\mathbf{r}, v) := -(\nabla \cdot \mathbf{r}, v)_h. \quad (5) $$

Notice that $c_h$ is continuous on $(\mathbf{P}_h \times V_h) \times (\mathbf{Q}_h \times V_h)$, for there is a norm $\|c_h\|$ of $c_h$ such that,

$$ c_h((\mathbf{r}; w), (\mathbf{q}; v)) \le \|c_h\| \|\| (\mathbf{r}; w) \|\|_h \|\| (\mathbf{q}; v) \|\|_h \quad \forall ((\mathbf{r}; w), (\mathbf{q}; v)) \in (\mathbf{P}_h \times V_h) \times (\mathbf{Q}_h \times V_h). \quad (6) $$
Actually we can take $\|c_h\| = \max[\|a_h\|, \|b_h\|, \|d_h\|]$ and observing that $\|a_h\| = \|b_h\| = \|d_h\| = 1$, we infer that $c_h$ is uniformly continuous, that is, its norm is independent of the particular subspaces $\mathbf{P}_h, \mathbf{Q}_h, V_h$ under consideration.

Now letting $F_h(v)$ be a suitable approximation of $(f, v)$ for any $v \in V_h$, the approximate problem to solve is

$$ \begin{cases} \text{Find } (\mathbf{p}_h; u_h) \in \mathbf{P}_h \times V_h \text{ such that,} \\ c_h((\mathbf{p}_h; u_h), (\mathbf{q}; v)) = F_h(v) \quad \forall (\mathbf{q}; v) \in \mathbf{Q}_h \times V_h. \end{cases} \tag{7} $$

Equation (7) is a standard linear variational problem defined upon a pair of different Hilbert spaces for the solution and the testing elements. Referring to [9], provided the right hand side functional is continuous, problems of this kind have a unique solution if and only if their bilinear form is continuous and weakly coercive over this pair of spaces. According to [16], in the case under study, weak coercivity means that

- There exists a constant $\gamma_h > 0$ such that $\forall (\mathbf{r}; w) \in \mathbf{P}_h \times V_h$ there is a non zero $(\mathbf{q}; v) \in \mathbf{Q}_h \times V_h$ such that $c_h((\mathbf{r}; w), (\mathbf{q}; v)) \ge \gamma_h |||(\mathbf{r}; w)|||_h |||(\mathbf{q}; v)|||_h$;

- The dimensions of $\mathbf{P}_h$ and $\mathbf{Q}_h$ are the same.

Now assume that there exists a function $\tilde{u}$ defined in a domain $\tilde{\Omega}$ large enough to contain both $\Omega$ and $\Omega_h$ for all meshes in use, which coincides with $u$ in $\Omega$ and fulfills $\tilde{u} \in H^1(\tilde{\Omega})$, $\Delta \tilde{u} \in L^2(\tilde{\Omega})$. The field $\tilde{\mathbf{p}} := \nabla \tilde{u}$ belongs to $\mathbf{H}(div, \tilde{\Omega})$.

Next, setting $\tilde{\Omega}_h := \Omega_h \cup \Omega \subset \tilde{\Omega}$, we observe that $c_h((\mathbf{r}; w), (\mathbf{q}; v))$ is well defined $\forall (\mathbf{r}; w) \in \mathbf{H}(div, \tilde{\Omega}_h) \times L^2(\tilde{\Omega}_h)$ and $\forall (\mathbf{q}; v) \in \mathbf{Q}_h \times V_h$ and thus $|||(\mathbf{r}; w)|||_h$ as well. Notice that $c_h$ is also continuous in the sense $c_h((\mathbf{r}; w), (\mathbf{q}; v)) \le \|c_h\| |||(\mathbf{r}; w)|||_h |||(\mathbf{q}; v)|||_h \forall ((\mathbf{r}; w); (\mathbf{q}; v)) \in [\mathbf{H}(div, \tilde{\Omega}_h) \times L^2(\tilde{\Omega}_h)] \times [\mathbf{Q}_h \times V_h]$ with $\|c_h\| = 1$.

On the other hand, in the case of curvilinear boundaries, in principle, we are dealing with an external approximation also called a variational crime (cf. [34]). Hence we resort to the general upper bound for linear variational problems given in [16] for weakly coercive bilinear forms defined upon two different spaces. In the case under study this writes

$$ \left\{ \begin{array}{l} |||(\tilde{\mathbf{p}}; \tilde{u}) - (\mathbf{p}_h; u_h)|||_h \le \frac{1}{\gamma_h} [\inf_{(\mathbf{r};w) \in \mathbf{P}_h \times V_h} |||(\tilde{\mathbf{p}}; \tilde{u}) - (\mathbf{r}; w)|||_h \\ + \sup_{(\mathbf{q};v) \in \mathbf{Q}_h \times V_h / |||(\mathbf{q};v)|||_h = 1} |c_h((\tilde{\mathbf{p}}; \tilde{u}), (\mathbf{q}; v)) - F_h(v)|]. \end{array} \right. \tag{8} $$

From (8) the constant $\gamma_h$ is seen to play a key role in error estimates for problem (7). Incidentally, the bilinear form $c_h$ is of the type studied in [30] with $P = \mathbf{P}_h$, $Q = \mathbf{Q}_h$ $U \equiv V = V_h$. Actually, as recalled in the introductory paragraph of Section 5, the key to establishing its weak coercivity is Theorem 3.1 of that article, whose assumptions are inf-sup conditions to be fulfilled by the bilinear forms $a_h, b_h$ and $d_h$ addressed in Sections 5.1, 5.2 and 5.3, respectively. Moreover it is proven in [30] that $\gamma_h$ can be expressed in terms of the constants of such conditions that supposedly hold for those bilinear forms. Naturally enough, the inf-sup conditions are established for the three of them in the aforementioned subsections, in connection with the particular form of (7) considered here.

More specifically, the spaces $\mathbf{Q}_h, \mathbf{P}_h$ and $V_h$ used in this work are linked to the Raviart-Thomas mixed finite element method of order $k$ for triangles [25]. We recall that (see e.g., [10]), if $\Omega$ is a polygon, provided $u$ belongs to $H^{k+2}(\Omega)$ (cf. [1]), this method is of order $k+1$ in the standard norm of the space $\mathbf{H}(div, \Omega) \times L^2(\Omega)$. Notice however that, in general, if $\Omega$ is a polygon, the above regularity can be expected at most if $k=0$ and $\Omega$ is convex. As a matter of fact, the purpose of problem (7) is to take the best advantage of higher order methods in this family, in the framework of problems posed in smooth domains, since in this case the exact solution can be much more regular.

### 3. MESHES AND RELATED SETS, NOTATIONS AND ASSUMPTIONS

The finite element spaces will shall work with are defined in connection with triangulations of $\Omega$, whose main features are next addressed.
<figure><img src="image_2.png" /><figcaption>FIGURE 1. Normal-flux DOFs shifted from $\Gamma_{1,h}$ to $\Gamma_1$ for $RT_1(T)$ (left) and $RT_0(T)$ (right) with $T \in S_{1,h}$.</figcaption></figure>

Let $\mathcal{F} := \{\mathcal{T}_h\}_h$ be a family of meshes consisting of straight-edged triangles, satisfying the well known compatibility conditions for the finite element method (cf. [13]), $h$ being the maximum edge length of all the triangles in the mesh $\mathcal{T}_h$. Every element of this mesh is to be viewed as a closed set and $\mathcal{T}_h$ is assumed to fit $\Omega$ in such a way that all the vertexes of the polygon $\Omega_h$ lie on $\Gamma$, where $\Omega_h$ is the interior of $\cup_{T \in \mathcal{T}_h} T$. We recall that $\tilde{\Omega}_h = \Omega \cup \Omega_h$ and further define $\Omega'_h := \Omega \cap \Omega_h$. The norms of $L^2(\Omega'_h)$, $\mathbf{H}(\mathrm{div}, \Omega'_h)$ and $\mathbf{H}(\mathrm{div}, \Omega'_h) \times L^2(\Omega'_h)$ are respectively denoted by $\| \cdot \|'_{0,h}$, $| | \cdot | |'_h$ and $||| \cdot ||'_h$.

We assume that the transition points between $\Gamma_1$ and $\Gamma_0$, if any, are always vertexes of the polygon $\Omega_h$. The boundary of $\Omega_h$ is denoted by $\Gamma_h$ and $n_h$ represents the unit outer normal vector to $\Gamma_h$. Finally we define $\Gamma_{1,h}$ to be the union of the edges of $\Gamma_h$ having both ends in the union of $\Gamma_1$ and the set of transition points between $\Gamma_0$ and $\Gamma_1$, if any, and let $\Gamma_{0,h}$ be the closure of $\Gamma_h \setminus \Gamma_{1,h}$.

$\mathcal{F}$ is assumed to be a regular family of partitions in the sense of [13] (cf. Sect. 3.1), though not necessarily quasi-uniform. The maximum edge length of every $T \in \mathcal{T}_h$ is denoted by $h_T$. Throughout this work we make the very reasonable assumption that no triangle $T \in \mathcal{T}_h$ has its three vertexes on $\Gamma$. It follows that every element in $\mathcal{T}_h$ has at most one edge contained in $\Gamma_h$.

Let $S_h$ (resp. $S_{1,h}$) be the subset of $\mathcal{T}_h$ consisting of triangles $T$ having one edge, say $e_T$, on $\Gamma_h$ (resp. $\Gamma_{1,h}$). We denote by $O_T$ the vertex of $T \in S_h$ not belonging to $\Gamma$ and by $n_T$ the unit outer normal vector along $e_T$. Notice that the interior of every triangle in $\mathcal{T}_h \setminus S_h$ has an empty intersection with $\Gamma_h$.

The following additional considerations are useful: $\forall T \in S_h$ we denote by $\Delta_T$ the closed subset of $\tilde{\Omega}_h$ not containing $T$ delimited by $\Gamma$ and the edge $e_T$ of $T$, as illustrated in Figure 1. Moreover, $\forall T \in S_h$, we define $T^\Gamma$ to be the set $T \cup \Delta_T$ and $T^{\Gamma'}$ to be the closure of the set $T \setminus \Delta_T$. Notice that, if $\Omega$ is convex $T^\Gamma = T$ and $\Omega_h$ is a proper subset of $\Omega$. Otherwise there are sets $T^\Gamma$ which coincide with $T$. Finally we set $\Gamma_T := T^\Gamma \cap \Gamma$ for any $T \in S_h$.

Although this is not a real restriction to the practical application of our method, for theoretical purposes and for the sake of simplicity, throughout this article we assume that every triangulation $\mathcal{T}_h \in \mathcal{F}$ is constructed in such a way that all the transition points between the convex and concave portions of $\Gamma$, if any, are vertexes of triangles thereof. In this respect, it is useful to define two subsets $S'_{0,h}$ and $S'_{1,h}$ of $S_{0,h}$ and $S_{1,h}$, which contain triangles having two vertexes on concave portions of $\bar{\Gamma}_0$ and $\bar{\Gamma}_1$, respectively. Notice that, owing to the above construction, such transition points are the same as those between convex and concave portions of $\Gamma_h$. We also set $S'_h := S'_{0,h} \cup S'_{1,h}$.

For convenience, the reliability analysis of our method will be carried out under the following assumption, which is by no means necessary, neither to define it, nor to implement it:

*Assumption*\*: $h$ is small enough for the intersection $N \in T^\Gamma$ with $\Gamma$ of the perpendicular to $e_T$ passing through any point $M \in e_T$ to be unique, $\forall T \in S_h$.
The main idea behind Assumption* is to ensure that $\Gamma_h$ follows very closely the shape of $\Gamma$, which might not be the case of a coarse mesh.

To close this section we prove:

**Proposition 3.1.** If Assumption* holds, there exists a constant $C_0$ depending only on $\Gamma$ such that $\forall M \in e_T$ the length of the segment perpendicular to $e_T$ joining $M$ and $N \in \Gamma_T$ is bounded above by $C_0h_T^2$.

For the proof of Propositions 3.1 we refer to Proposition 2.1 of [27] or yet to [23], Lemma A.2.

**Proposition 3.2.** If $\Gamma$ is of the piecewise $C^2$-class and Assumption* holds, $\forall T \in S_{1,h}$ and $\forall N \in \Gamma_T$ the outer normal $\mathbf{n}(N)$ to $\Gamma_T$ at $N$ satisfies $|\mathbf{n}(N) - \mathbf{n}_T| \le C_1h_T$, where $C_1$ is a constant depending only on $\Gamma_1$.

*Proof.* First we recall Proposition 3.2 of [27], according to which there is a constant $C_\Gamma$ depending only on $\Gamma$ such that

$$|\tan \varphi| \le C_\Gamma h_T, \text{ where } \varphi \text{ is the angle between } \mathbf{n}(N) \text{ and } \mathbf{n}_T. \quad (9)$$

On the other hand we have $|\mathbf{n}(N) - \mathbf{n}_T|^2 = 2(1 - \mathbf{n}(N) \cdot \mathbf{n}_T) = 2(1 - \cos\varphi)$. Thus $|\mathbf{n}(N) - \mathbf{n}_T|^2 = 2\sin^2\varphi/(1 + \cos\varphi) = 2\tan^2\varphi/(\cos^{-2}\varphi + \cos^{-1}\varphi)$, that is, $|\mathbf{n}(N) - \mathbf{n}_T|^2 = 2\tan^2\varphi/(\tan^2\varphi + 1 + \sqrt{1 + \tan^2\varphi}) \le \tan^2\varphi$. It trivially follows that $|\mathbf{n}(N) - \mathbf{n}_T| \le C_\Gamma h_T$, i.e., $C_1 = C_\Gamma$. $\square$

## 4. THE FINITE-ELEMENT APPROXIMATION METHOD

In this section we give the specific form of the finite element approximation spaces $V_h$, $Q_h$ and $P_h$ in (7), associated with $\mathcal{T}_h$ and an integer $k \ge 0$, that we advocate for the numerical solution of (1).

First of all we denote by $\text{RT}_k(T)$ the Raviart–Thomas space of vector fields defined in a triangle $T$ (cf. [25]), whose dimension is $(k+1)(k+3)$.

- $V_h$ is the space of functions $v$ such that $\forall T \in \mathcal{T}_h, v|_T \in P_k(T)$, satisfying $v \in L_0^2(\Omega_h)$ if $\text{length}(\Gamma_0) = 0$.

- $Q_h$ is the space consisting of fields $q \in \mathbf{H}(\text{div}, \Omega_h)$, whose restriction to each triangle $T \in \mathcal{T}_h$ belongs to $\text{RT}_k(T)$, and whose normal component vanishes on $\Gamma_{1,h}$. We recall that, classically, the DOFs of $\text{RT}_k(T)$ are the $k+1$ moments of order $k$ of the field's normal components along each one of the edges of $T$ and, for $k > 0$ only, the $k(k+1)$ moments of order $k-1$ in $T$ of both field's components. For convenience we extend every field $q \in Q_h$ by zero in $\bar{\Omega} \setminus \bar{\Omega}_h$.

- $P_h$ in turn is the space of fields $p$ defined in $\tilde{\Omega}_h$ having the properties listed below.

(1) The restriction of $p$ to every $T \in \mathcal{T}_h$ belongs to the Raviart–Thomas space $\text{RT}_k(T)$;

(2) The normal component of every $p$ is continuous across all the edges of the elements in $\mathcal{T}_h$ contained in the interior of $\Omega_h$;

(3) In $\bar{\Omega} \setminus \bar{\Omega}_h$ a field $p$ is expressed in such a way that $\forall T \in S_h$ its polynomial expression in $T$ extends to $\Delta_T$ in case $\Delta_T \setminus T$ is not empty;

(4) $\forall T \in S_{1,h}$, the (outer) normal component of $p$ vanishes at every point which is the nearest intersection with $\Gamma_1$ of the perpendicular to the edge $e_T$ passing through one of the $k+1$ Gauss-quadrature points of $e_T$. The intersection point on $\Gamma_T$ and the corresponding Gauss-quadrature point are generically denoted by $N$ and $M$, respectively, as illustrated in Figure 1.

Before establishing that $P_h$ is a non empty finite-dimensional space we recall the following lemma, whose proof is given in [27], except for the obvious left inequality in (11):

**Lemma 4.1** ([27]). Provided $h$ satisfies Assumption*, there exist three constants $C_\infty$, $C_0^-$ and $C_0^+$ depending only on $\Gamma$ and the shape regularity of $\mathcal{T}_h$ (cf. [8], Ch.4, Sect. 4) such that $\forall w \in P_k(T^\Gamma)$ and $\forall T \in S_h$ it holds:

$$\|w\|_{0,\infty,T^\Gamma} \le C_\infty \|w\|_{0,\infty,T^{\Gamma'}} \quad (10)$$

and

$$C_0^- h_T^{-1} \|w\|_{0,T^\Gamma} \le \|w\|_{0,\infty,T^\Gamma} \le C_0^+ h_T^{-1} \|w\|_{0,T^{\Gamma'}}. \quad (11)$$
Now we state

**Proposition 4.2.** Provided $h$ is small enough, given an arbitrary set $B_k$ of $m_k$ real values with $m_k = (k + 1)(k + 2), \forall T \in S_{1,h}$ there exists a unique field in $RT_k(T)$, whose normal components across $\Gamma_1$ vanish at the $k + 1$ points $N \in \Gamma_T$ defined in accordance with item 4. of the above definition of $P_h$, and whose remaining $m_k$ DOFs of fields in the space $RT_k(T)$ take up their values in the set $B_k$ on a one-to-one basis.

*Proof.* Let $B_k = \{b_1, b_2, \dots, b_{m_k}\}$ and $b$ be the vector of $\Re^{n_k}$ given by $b := [0, \dots, 0, b_1, \dots, b_{m_k}]$, i.e., $b$ is a vector whose $k + 1$ first entries are zero. If the boundary nodes $N$ were replaced by the corresponding $M \in \Gamma_{1,h} \cap T$, it is clear that the result would hold true, according to the well known properties of the space $RT_k(T)$. The vector $a$ of coefficients $a_i$ for $i = 1, 2, \dots, n_k$ of the canonical basis fields $q_{i,T}$ of $RT_k(T)$ for $1 \le i \le n_k$ would be given by $a_i = b_i$ for $1 \le i \le n_k$.

Let us denote by $\mathcal{H}_i$ the DOFs of $RT_k(T)$ and by $q_{i,T}$ the underlying canonical basis fields, for $i = 1, 2, \dots, n_k$, in such a way that the first $k + 1$ ones correspond to normal components across $e_T$. In doing so, the matrix $E$ whose entries are $e_{ij} := \mathcal{H}_i(q_{j,T})$ is the identity matrix. Now we define $\tilde{\mathcal{H}}_i$ to be $\mathcal{H}_i$ if $k + 1 < i \le n_k$ and otherwise to be the normal component across $\Gamma_1$ at one of the $k + 1$ points $N$ constructed in accordance with item 4. of the definition of $P_h$. For convenience we further number these points $N$ and the corresponding point $M$ on $e_T$, respectively by $N_i$ and $M_i$, for $i = 1, \dots, k + 1$. We also denote by $n_i$ the outer normal to $\Gamma_1$ at $N_i$, for $i = 1, 2, \dots, k + 1$.

The proposition will be proved if the $n_k \times n_k$ linear system $\tilde{E}a = b$ is uniquely solvable, where $\tilde{E}$ is the matrix with entries $\tilde{e}_{ij} := \tilde{\mathcal{H}}_i(q_{j,T})$. Clearly we have $\tilde{E} = E + D_E$, where the entries of $D_E$ are $d_{ij} := \tilde{\mathcal{H}}_i(q_{j,T}) - \mathcal{H}_i(q_{j,T})$. While $d_{ij} = 0 \forall j$ if $k + 1 < i \le n_k$, in case $1 \le i \le k + 1$ we have

$$d_{i,j} = q_{j,T}(N_i) \cdot n_i - q_{j,T}(M_i) \cdot n_T.$$

This means that,

$$ \left\{ \begin{array}{l} d_{ij} = r_{ij} + s_{ij}, \text{ where} \\ r_{ij} = q_{j,T}(N_i) \cdot n_i - q_{j,T}(N_i) \cdot n_T \\ s_{ij} = q_{j,T}(N_i) \cdot n_T - q_{j,T}(M_i) \cdot n_T. \end{array} \right. \qquad (12) $$

From Proposition 3.2 we may write $|n_i - n_T| \le C_T h_T$ for $i = 1, 2, \dots, k + 1$. Taking into account the fact that there is a mesh-independent constant $\tilde{C}_q$ such that $\|q_{j,T}\|_{0,\infty,T} \le \tilde{C}_q$ for $i = 1, 2, \dots, k + 1$, it follows that,

$$|r_{ij}| \le C_R h_T \quad \forall i, j, \text{ with } C_R = \tilde{C}_q C_T. \qquad (13)$$

On the other hand, from Proposition 3.1, the length of the segment $\overline{M_i N_i}$ is bounded above by $C_0 h_T^2$. Thus, a simple Taylor expansion of the function $q_{j,T} \cdot n_T$ implies that,

$$|s_{ij}| \le C_0 h_T^2 \max_{x \in T^\Gamma} |\nabla\{q_{j,T} \cdot n_T\}(x)|,$$

or yet, recalling (10),

$$|s_{ij}| \le C_0 C_\infty h_T^2 \max_{x \in T} |\nabla\{q_{j,T} \cdot n_T\}(x)|. \qquad (14)$$

Finally, from Lemma 4.1 and a standard argument we know that the maximum in (14) is bounded above by a mesh-independent constant times $h_T^{-1}$. In short we have

$$|s_{ij}| \le C_S h_T \quad \forall i, j, \qquad (15)$$

where $C_S$ is a constant independent of $T$.

It follows from (12), (13) and (15) that the matrix $\tilde{E}$ equals the identity matrix plus an $O(h_T)$ matrix $D_E = [d_{ij}]$. Therefore $\tilde{E}$ is an invertible matrix, as long as $h$ is sufficiently small, more precisely, if $h < [\tilde{C}_k]^{-1/2}$, with $\tilde{C}_k = (C_R + C_S)/\Theta$ for a suitably small $\Theta$. $\square$
Thanks to Proposition 4.2 it is legitimate to pose problem (7) for the spaces $\mathbf{P}_h$, $\mathbf{Q}_h$ and $V_h$ defined in Section 3.1, whose solution $(\mathbf{p}_h; u_h)$ approximates the solution of (1) with $\dim \mathbf{P}_h = \dim \mathbf{Q}_h$.

We close this section by proving a property of space $\mathbf{P}_h$ to be repeatedly used in the sequel.

**Lemma 4.3.** Assume that $\Gamma$ is of the $C^{2k+2}$-class and of the piecewise $C^{2k+3}$-class.

Given $\mathbf{p}$ in $\mathbf{P}_h$ let $g$ be its normal trace on $\Gamma_1$. Let $v$ be given in $V_h$. There exists a constant $C_{\Gamma_1}$ such that $\forall T \in S_{1,h}$ we have

$$ \left| \int_{\Gamma_T} g v \, ds \right| \le C_{\Gamma_1} h_T^{k+1} \| \mathbf{p} \|_{k+1,T} \| v \|_{0,T}. \quad (16) $$

*Proof.* First of all we denote by $\phi$ the function of $x \in [0, l_T]$ representing the $y$-coordinate of a point $N \in \Gamma_T$ corresponding to $M \in e_T$ with abscissa $x$ in the local direct cartesian coordinate system $(A_T, x, y)$, whose origin $A_T$ is chosen in such a manner that $y(x) \ge 0 \ \forall x$. In doing so, for any bounded function $a$ defined in a neighborhood of $\Gamma_T$, we define $\check{a}$ in $[0, l_T]$ by $\check{a}(x) = a[x, \phi(x)]$. It is clear that there exists a mesh-independent constant $C_q$ such that for every such functions $a$ we have

$$ \left| \int_{\Gamma_T} a \, ds \right| \le C_q \left| \int_0^{l_T} \check{a} \, dx \right|. \quad (17) $$

According to our regularity assumption on $\Gamma$ the derivative of order $2k+3$ of $\phi$ and the derivative of order $2k+2$ of the normal vector $\mathbf{n}$ are uniformly bounded in $(0, l_T)$. Hence $\check{g}$ belongs to $W^{2k+2,\infty}(0, l_T)$.

Now we consider $a$ to be the function $gv$. Noticing that, by construction, $\check{g}$ vanishes at the $k+1$ Gauss-quadrature points of $e_T$, so does $\check{a}$. Thus by well known properties of Gauss integration in one-dimensional space (see e.g., [7, 24]), it holds for a constant $\check{C}$ independent of $T$,

$$ \left| \int_{0}^{l_T} \check{a} \, dx \right| \le \check{C} l_T^{2k+3} \sup_{x \in (0, l_T)} \left| \frac{d^{2k+2} \check{a}}{dx^{2k+2}}(x) \right|. \quad (18) $$

Combining (17) and (18), setting $C_a := C_q \check{C}$ we easily come up with

$$ \left| \int_{\Gamma_T} a \, ds \right| \le C_a h_T^{2k+3} |\check{a}|_{2k+2,\infty, e_T}. \quad (19) $$

Since $\check{a} = \check{g} \check{v}$, after straightforward manipulations we obtain for another mesh-independent constant $C'_a$

$$ \left| \int_{\Gamma_T} a \, ds \right| \le C'_a h_T^{2k+3} \|\check{g}\|_{2k+2,\infty, e_T} \|\check{v}\|_{2k+2,\infty, e_T}. \quad (20) $$

On the other hand, we note that $\mathbf{n}$ can be expressed only in terms of $x$ in a neighborhood of $\Gamma_T$ and moreover $\mathbf{n}_{|(A_T,x,y)} = (-\sin \varphi; \cos \varphi)$. Hence we can uniformly bound above the derivatives of $\mathbf{n}$ up to the $2k+2$-th order as follows.

Denoting by $\mathcal{A}^{(j)}(x)$ the $j$-th order derivative for $j > 3$ of a function or vector field $\mathcal{A}(x)$, we have

$$ \mathbf{n}' = (-\cos \varphi; -\sin \varphi)\varphi' $$

where $\varphi$ is defined by $\tan \varphi := \phi'$, with $\varphi \in [-\pi/2, \pi/2]$. Therefore

$$ \varphi' = \frac{\phi''}{[1 + (\phi')^2]^2}. $$

Noticing that $|\phi'|$ is uniformly bounded above by a constant times $h_T$ in $e_T$ and the functions $|\phi^{(j)}|$ are all uniformly bounded above by mesh-independent constant $C_\phi^{(j)}$ for $2 \le j \le 2k+2$, we easily obtain

$$ |\mathbf{n}'| \le C_{\mathbf{n}}^{(1)} := C_{\phi}^{(2)}. \quad (21) $$
Performing repeatedly the same type of calculations for higher order derivatives of $\mathbf{n}$, we conclude with no difficulty that the moduli of the derivatives $\mathbf{n}^{(j)}$ of any order up to the $2k+2$-th are uniformly bounded above independently of $T$ by constants $C_{\mathbf{n}}^{(j)}$, expressed in terms of the constants $C_{\phi}^{(i)}$ for $i = 1, 2, \dots, j+1$ that is

$$|\mathbf{n}^{(j)}| \le C_{\mathbf{n}}^{(j)} \quad \text{for } j = 1, 2, \dots, 2k+2. \qquad (22)$$

Now using the chain rule, we have

$$\check{g}' = \left[ \frac{\partial \mathbf{p}}{\partial x} + \frac{\partial \mathbf{p}}{\partial y} \phi' \right] \cdot \mathbf{n} + \mathbf{p} \cdot \mathbf{n}'. \qquad (23)$$

Then the following bound readily derives from (21) and (23):

$$|\check{g}'(x)| \le C_g^{(1)} \sqrt{|[\nabla \mathbf{p}](x)|^2 + |\mathbf{p}(x)|^2} \quad \forall x \in [0, l_T] \qquad (24)$$

with $C_g^{(1)} := \sqrt{1 + [C_{\phi}^{(1)}]^2 + [C_{\mathbf{n}}^{(1)}]^2}$.

Further differentiating $\check{g}$, taking into account (22) and noticing that $D^i \mathbf{p}$ vanishes identically for $i > k + 1$, akin to Theorem 4.4 of [27] we find out that there exists a constant $C_g^{(2k+2)}$ independent of $T$ expressed in terms of $C_{\mathbf{n}}^{(j)}$ and $C_{\phi}^{(j)}$ such that, for $j = 1, 2, \dots, 2k+2$,

$$|\check{g}^{(j)}(x)| \le C_g^{(2k+2)} \left\{ \sum_{i=0}^{k+1} |[D^i \mathbf{p}](x)|^2 \right\}^{1/2} \quad \forall x \in [0, l_T]. \qquad (25)$$

As for $\check{v}$, in the light of the above calculations for $\check{g}$, and observing that $D^i v$ vanishes identically for $i > k$, it is even easier to infer the existence of a mesh-independent constant $C_v^{(2k+2)}$ expressed in terms of $C_{\phi}^{(j)}$ such that, for $j = 1, 2, \dots, 2k+2$,

$$|\check{v}^{(j)}(x)| \le C_v^{(2k+2)} \left\{ \sum_{i=0}^{k} |[D^i v](x)|^2 \right\}^{1/2} \quad \forall x \in [0, l_T]. \qquad (26)$$

Plugging (25) and (26) into (20) we obtain,

$$\int_{\Gamma_T} \mathbf{p} \cdot \mathbf{n} \, v \, ds \le C'_a C_g^{(2k+2)} C_v^{(2k+2)} (k+1) h_T^{2k+3} \|\mathbf{p}\|_{k+1, \infty, T^\Gamma} \|v\|_{k, \infty, T^\Gamma}. \qquad (27)$$

Applying (11) to each norm in (27) and using $k$ times a classical inverse inequality (see e.g., [19]) to bound above $\|v\|_{k,T}$ by $h_T^{-k} \|v\|_{0,T}$ multiplied by a constant, we immediately come up with (16) for a suitable constant $C_{\Gamma_1}$. $\square$

## 5. WELL-POSEDNESS AND UNIFORM STABILITY OF THE APPROXIMATE PROBLEM

Besides establishing that problem (7) has a unique solution, as a consequence of the fact that the bilinear form $c_h$ given by (4) and (5) is weakly coercive, we will prove in Theorems 5.8 and 5.13 hereafter that its weak coercivity constant $\gamma_h$ is mesh-independent. This implies the uniform stability of (7). According to [30], this result is a consequence of the validity of uniform inf-sup conditions for the bilinear forms $a_h$, $b_h$ and $d_h$ defined in (5). Each one of them is worked out in one of the three subsections that follow.

At the end of Section 5.3 the inf-sup constant $\gamma_h$ of $c_h$ will be denoted by either $\gamma$ or $\gamma'$, depending on whether $\Omega$ is convex or not. Both are exhibited therein as uniform stability constants derived from the resulting uniform inf-sup constants of $a_h$, $b_h$ and $d_h$ denoted respectively by $\alpha$, $\beta$ and $\delta$ or $\delta'$ accordingly. More precisely the corresponding inf-sup conditions are (29), (42) and (87) (or yet (142)) stated in Theorems 5.1, 5.3 and 5.7 (or yet 5.12).
**5.1. Uniform inf-sup condition for bilinear form $a_h$ on $Ker(div)$**

**Theorem 5.1.** Let

$$
\begin{align}
\mathbf{R}_h &:= \{\mathbf{r} | \mathbf{r} \in \mathbf{P}_h, \text{ and } d_h(\mathbf{r}, v) = 0 \quad \forall v \in V_h\} \text{ and} \tag{28} \\
\mathbf{S}_h &:= \{\mathbf{s} | \mathbf{s} \in \mathbf{Q}_h, \text{ and } b_h(v, \mathbf{s}) = 0 \quad \forall v \in V_h\}. \nonumber
\end{align}
$$

Provided $h$ is sufficiently small, there exists a constant $\alpha > 0$ independent of $h$ such that,

$$
\forall \mathbf{r} \in \mathbf{R}_h \exists \mathbf{s} \in \mathbf{S}_h \setminus \{\mathbf{0}\} \text{ such that } \frac{a_h(\mathbf{r}, \mathbf{s})}{[|\mathbf{s}|]_h} \ge \alpha [|\mathbf{r}|]_h. \quad (29)
$$

Proof. It is easy to see that the divergence of any field in RT$_k(T)$ is a polynomial in $P_k(T)$ whose $k+1$ terms of degree $k$ have the same coefficients as those of the terms of degree $k+1$ of the field itself multiplied by $k+2$. As a matter of fact, the divergence operator is a surjection from RT$_k(T)$ onto $P_k(T)$. Hence the divergence of every field in either $\mathbf{R}_h$ or $\mathbf{S}_h$ vanishes identically. We also observe that the restriction to $T$ of a field in either $\mathbf{R}_h$ or $\mathbf{S}_h$ is a field in a subspace RT'$_k(T)$ of $[P_k(T)]^2$, whose dimension is $2dimP_k(T) - k(k+1)/2$, that is, $l_k := (k+4)(k+1)/2$.

For $T \in \mathcal{T}_h$ let $\{\mathbf{q}'_{i,T}\}_{i=1}^{l_k}$ be a basis of RT'$_k(T)$ associated with $l_k$ DOFs $\{\mathcal{G}_{i,T}(\mathbf{q}'_{j,T})\}_{i=1}^{l_k}$ such that $\mathcal{G}_{i,T}(\mathbf{q}'_{j,T}) = \delta_{ij}$.
Since the normal component of a field in $\mathbf{S}_h$ vanishes identically on $\Gamma_{1,h}$, for every $T \in \mathcal{S}_{1,h}$ we choose these
DOFs in such a way that the first $k+1$ ones are the field's outer normal components along the edge $e_T$ at the
$k+1$ Gauss-quadrature points $M_i$ of $e_T$.

Now given a field $\mathbf{r} \in \mathbf{R}_h$, let $\mathbf{s} \in \mathbf{S}_h$ be such that

$$
- \forall T \in \mathcal{T}_h \setminus \mathcal{S}_{1,h} \quad \mathbf{s}_{|T} \equiv \mathbf{r}_{|T};
$$

$$
- \forall T \in \mathcal{S}_{1,h} \mathcal{G}_{i,T}(\mathbf{s}_{|T}) = \mathcal{G}_{i,T}(\mathbf{r}_{|T}), i = k+2, \dots, l_k.
$$

Since $\mathbf{s}$ only differs from $\mathbf{r}$ in a triangle $T \in \mathcal{S}_{1,h}$, we have

$$
\mathbf{s}_{|T} - \mathbf{r}_{|T} = \sum_{i=1}^{k+1} [\mathcal{G}_{i,T}(\mathbf{s}_{|T}) - \mathcal{G}_{i,T}(\mathbf{r}_{|T})] \mathbf{q}'_{i,T}. \quad (30)
$$

Abusively denoting $\mathbf{r}_{|T\mathrm{r}}$ by $\mathbf{r}_{|T}$, we note that $\mathcal{G}_{i,T}(\mathbf{s}_{|T}) = 0$ and that $\mathcal{G}_{i,T}(\mathbf{r}_{|T}) = \mathbf{r}_{|T}(M_i) \cdot \mathbf{n}_{T}$, for $1 \le i \le k1$.
In view of this and taking into account that $\mathbf{r}_{|T}(N_i) \cdot \mathbf{n}_i = 0$, we have,

$$
\mathcal{G}_{i,T}(\mathbf{s}_{|T}) - \mathcal{G}_{i,T}(\mathbf{r}_{|T}) = [\mathbf{r}_{|T}(M_i) \cdot \mathbf{n}_T - \mathbf{r}_{|T}(N_i) \cdot \mathbf{n}_T] + [\mathbf{r}_{|T}(N_i) \cdot \mathbf{n}_T - \mathbf{r}_{|T}(N_i) \cdot \mathbf{n}_i]. \quad (31)
$$

Now we resort to Propositions 3.1 and 3.2 to write,

$$
|\mathcal{G}_{i,T}(\mathbf{s}_{|T}) - \mathcal{G}_{i,T}(\mathbf{r}_{|T})| \le |\mathbf{r}_{|T}(M_i) - \mathbf{r}_{|T}(N_i)| + |\mathbf{r}_{|T}(N_i)||\mathbf{n}_T - \mathbf{n}_i| \le C_0 ||\nabla \mathbf{r}||_{0,\infty,T\mathrm{r}} h_T^2 + C_\Gamma h_T ||\mathbf{r}||_{0,\infty,T\mathrm{r}}. \quad (32)
$$

From Lemma 4.1 and a classical inverse inequality, we infer the existence of a mesh-independent constant $C_G$
such that

$$
|\mathcal{G}_{i,T}(\mathbf{s}_{|T}) - \mathcal{G}_{i,T}(\mathbf{r}_{|T})| \le C_G \|\mathbf{r}_{|T}\|_{0,T}. \quad (33)
$$

On the other hand, taking into account that there exists a mesh-independent constant $C'_{\mathbf{q}}$ such that $\|\mathbf{q}'_{i,T}\|_{0,T} \le C'_{\mathbf{q}}h_T$ for $i=1, \dots, k+1$, combining (30) and (33) we easily obtain with $C'_{k} := (k+1)C_G C'_{\mathbf{q}}$:

$$
\| \mathbf{s}_{|T} - \mathbf{r}_{|T} \|_{0,T} \le C'_{k} \| \mathbf{r}_{|T} \|_{0,T} h_T, \quad (34)
$$

which immediately yields:

$$
[|\mathbf{s} - \mathbf{r}|]_h \le C'_k h [|\mathbf{r}|]_h. \quad (35)
$$

From (35) we easily conclude that $[|\mathbf{s}|]_h \le (1 + C'_k h)[|\mathbf{r}|]_h$ and $[\mathbf{s}, \mathbf{r}]_h \ge (1 - C'_k h)[|\mathbf{r}|]_h^2$. Hence, recalling Proposition 4.2, for $h \le \max[\hat{C}_k, C'_k]^{-1/2}$, the result follows with $\alpha = 1/3$. $\square$

In the sequel we will use (29) in the following equivalent form:

$$
\forall \mathbf{r} \in \mathbf{R}_h \exists \mathbf{s} \in \mathbf{S}_h \text{ with } [|\mathbf{s}|]_h = 1 \text{ such that } a_h(\mathbf{r}, \mathbf{s}) \ge \alpha [|\mathbf{r}|]_h. \quad (36)
$$
## 5.2. Uniform inf-sup condition for bilinear form $b_h$

In order to prove this second crucial result for establishing the uniform stability of problem (7), namely, conditions (42) hereafter, we need the following preparatory material:

Similarly to [10], we make use of the $L^2(\Omega_h)$-projection operator $\pi_h$ onto the space $V_h$ and consider as well the interpolation operator $\Pi_h$ onto the space $Q_h$. The latter applies to fields in the space $W_h$ defined in connection with a real number $s$ strictly greater than two, i.e.,

$$W_h := \{\mathbf{q} | \mathbf{q} \in \mathbf{H}(\mathrm{div}, \Omega_h) \cap [L^s(\Omega_h)]^2, \mathbf{q} \cdot \mathbf{n}_h = 0 \text{ on } \Gamma_{1,h}\},$$

$W_h$ is equipped with the norm $\|\cdot\|_{W_h} := \|\cdot\|_{0,s,h} + \|\nabla \cdot\|_{0,h}$.

Let $\bar{V}$ be the space $L^2(\Omega_h)$ if $\mathrm{length}(\Gamma_0) > 0$ or $L_0^2(\Omega_h)$ otherwise.

**Proposition 5.2.** There exists a constant $\beta_W > 0$ independent of $h$ such that,

$$\forall v \in \bar{V} \quad \sup_{\mathbf{q} \in W_h} b_h(v, \mathbf{q}) \ge \beta_W \|v\|_{0,h}. \quad (37)$$

*Proof.* Let $v$ be given in $\bar{V}$ and $\mathbf{g} := -\nabla z$, where $z$ is the unique solution of the Poisson problem

$$\left\{
\begin{array}{ll}
-\Delta z = v & \text{in } \Omega_h \\
z = 0 & \text{on } \Gamma_{0,h} \text{ if } \mathrm{length}(\Gamma_{0,h}) \neq 0 \text{ and } \int_{\Omega_h} z = 0 \text{ otherwise} \\
\partial z / \partial n_h = 0 & \text{on } \Gamma_{1,h}.
\end{array}
\right.
\quad (38)$$

Even if $\mathrm{length}(\Gamma_0) > 0$, according to [22], $z \in H^r(\Omega_h)$ for some $r \in (3/2, 2)$ and moreover by the Rellich-Kondrachev Theorem [1], $\nabla z \in L^s(\Omega_h)$ for some $s > 2$. Actually, due to both arguments we know that there exists a constant $C_{s,h}$ depending on $\Omega_h$ such that

$$\|\nabla z\|_{0,s,h} \le C_{s,h} \|v\|_{0,h}. \quad (39)$$

It trivially follows that

$$\frac{b_h(v, \mathbf{g})}{\|\mathbf{g}\|_{0,s,h}} \ge C_{s,h}^{-1} \|v\|_{0,h} \quad (40)$$

and, as one can easily check,

$$\frac{b_h(v, \mathbf{g})}{\|\mathbf{g}\|_{W_h}} \ge \beta_h \|v\|_{0,h}, \quad (41)$$

with $\beta_h = (C_{s,h} + 1)^{-1}$.

In Appendix I of [4] we claim and prove that the constant $C_{s,h}$ has an upper bound $C_s$ independent of $h$. Hence the constant $\beta_h$ has a lower bound $\beta_W > 0$ independent of $h$. Setting $\tilde{\mathbf{q}} = \mathbf{g}/\|\mathbf{g}\|_{W_h}$ (37) follows. $\square$

Now we are able to establish the uniform inf-sup condition for $b_h$ in

**Theorem 5.3.** There exists a constant $\beta > 0$ independent of $h$ such that

$$\forall v \in V_h \quad \sup_{\mathbf{q} \in Q_h} b_h(v, \mathbf{q}) \ge \beta \|v\|_{0,h}. \quad (42)$$

*Proof.* According to (37) we know that the divergence operator has a continuous lifting from $\bar{V}$ to $W_h$. Moreover, following the same steps as in Chapter III of [10], we can assert that there exists a constant $C_\Pi$ independent of $h$ such that

$$[\Pi_h \mathbf{q}]_h \le C_\Pi \|\mathbf{q}\|_{W_h} \quad \forall \mathbf{q} \in W_h.$$
On the other hand, quoting again [10], for every $v \in V_h$ there exists $\mathbf{q} \in \mathbf{Q}_h$ such that $\nabla \cdot \mathbf{q} = v$. After a quick check, all the other conditions specified on page 138, Chapter IV of [10] for the commuting diagram (I.23) given therein, apply to the pairs of spaces $\mathbf{W}_h - \bar{V}$ on the upper part and $\mathbf{Q}_h - V_h$ on the lower part, while $\Pi_h$ and $\pi_h$ operate from one space to the other on its left and right parts.

Summarizing, $\forall \mathbf{q} \in \mathbf{W}_h$ it holds:

$$ \left\{ \begin{array}{l} b_h(v, \Pi_h \mathbf{q}) = b_h(v, \mathbf{q}) \quad \forall v \in V_h \\ [\Pi_h \mathbf{q}]_h \le C_{\Pi} \|\mathbf{q}\|_{\mathbf{W}_h}. \end{array} \right. \qquad (43) $$

Now, taking an arbitrary $v \in V_h$ and recalling (37), let $\mathbf{q} \in \mathbf{W}_h$ be such that $\|\mathbf{q}\|_{\mathbf{W}_h} = 1$ and $b_h(v, \mathbf{q}) \ge \beta \mathbf{w} \|v\|_{0,h}$. Then, from (43) we immediately conclude that (42) holds with $\beta = \beta_{\mathbf{w}}/C_{\Pi}$. $\square$

### 5.3. Uniform inf-sup conditions for bilinear forms $d_h$ and $c_h$

Taking into account the results of Sections 4.1 and 4.2, according to [30] a uniform inf-sup holds for bilinear form $c_h$, as long as a uniform inf-sup condition also holds for bilinear form $d_h$. In order to establish the latter, the following additional preparatory material is of paramount importance.

First of all we introduce for convenience the subset $\mathcal{F}'$ of triangulations of interest in the family $\mathcal{F}$. Henceforth the latter expression will qualify meshes $\mathcal{T}_h$ which are not too coarse. We refrain from further elaborating about this concept, since it implicitly goes hand in hand with an argument repeatedly invoked throughout this paper, namely, the assumption that $h$ is sufficiently small.

We shall need an interpolation operator onto the space $\mathbf{P}_h$ for fields defined beyond $\Omega_h$. The other way around, fields defined in $\Omega$ are not necessarily defined in $\Omega_h$. Thus, in the remainder of this section and on we will work with extensions thereof outside $\Omega$, say, in a smooth domain $\tilde{\Omega}$ strictly containing $\tilde{\Omega}$. More specifically we shall use the following auxiliary material:

We require $\tilde{\Omega}$ to be sufficiently large to strictly contain the sets $\tilde{\Omega}_h$ for all $\mathcal{T}_h \in \mathcal{F}'$. We denote by $\tilde{\Gamma}$ the boundary of $\tilde{\Omega}$ and define $\Delta_{\tilde{\Omega}} = \tilde{\Omega} \setminus \tilde{\Omega}$. Notice that the boundary $\partial\Delta_{\tilde{\Omega}}$ of $\Delta_{\tilde{\Omega}}$ equals $\Gamma \cup \tilde{\Gamma}$.

Now for any $\mathbf{m} \in [W^{t,r}(\tilde{\Omega})]^2 \cap \mathbf{H}(\mathrm{div}, \tilde{\Omega})$, with $t \in \mathbb{R}^+$ and $r \ge 1$, let us assume that there exists an extension $\tilde{\mathbf{m}}$ of $\mathbf{m}$ to the whole $\tilde{\Omega}$ such that $\tilde{\mathbf{m}} \in [W^{t,r}(\tilde{\Omega})]^2 \cap \mathbf{H}(\mathrm{div}, \tilde{\Omega})$. Referring to an argument highlighted in [10], as long as $L^s(\tilde{\Omega})$ is continuously embedded into $W^{t,r}(\tilde{\Omega})$ for some $s > 2$, we may define an interpolation operator $\tilde{\Pi}_h : [W^{t,r}(\tilde{\Omega})]^2 \cap \mathbf{H}(\mathrm{div}, \tilde{\Omega}) \to \mathbf{P}_h$ by the following relations, where $\partial T$ represents the boundary of a triangle $T \in \mathcal{T}_h$:

$$ \left\{ \begin{array}{ll} \forall T \in \mathcal{T}_h \setminus \mathcal{S}_{1,h}, \quad \int_e [\tilde{\Pi}_h \tilde{\mathbf{m}} - \tilde{\mathbf{m}}] \cdot \mathbf{n}_T v = 0 & \forall v \in P_k(e) \quad \forall \text{ edge } e \subset \partial T; \\ \forall T \in \mathcal{S}_{1,h} \quad \int_e [\tilde{\Pi}_h \tilde{\mathbf{m}} - \tilde{\mathbf{m}}] \cdot \mathbf{n}_T v = 0 & \forall v \in P_k(e) \quad \forall \text{ edge } e \subset \partial T \setminus e_T; \\ \text{For } k > 0, \quad \int_T [\tilde{\Pi}_h \tilde{\mathbf{m}} - \tilde{\mathbf{m}}] \cdot \mathbf{o} = 0 & \forall \mathbf{o} \in [P_{k-1}(T)]^2. \end{array} \right. \qquad (44) $$

We can assert that $\tilde{\Pi}_h \tilde{\mathbf{m}} \in \mathbf{H}(\mathrm{div}, \tilde{\Omega}) \forall \tilde{\mathbf{m}} \in [W^{t,r}(\tilde{\Omega})]^2 \cap \mathbf{H}(\mathrm{div}, \tilde{\Omega})$. Moreover, as long as for every bounded Lipschitz domain $D$ of $\mathbb{R}^2$, $W^{t,r}(D)$ is continuously embedded into $W^{q,s}(D)$ for some $q \ge 0$ and $s \ge 1$, there exists a constant $C_{q,s}^{t,r}$ not depending on $T$ such that the following estimate holds

$$ \| \tilde{\Pi}_h \tilde{\mathbf{m}} - \tilde{\mathbf{m}} \|_{q,s,T^r} \le C_{q,s}^{t,r} h_T^{2/s-2/r} h_T^{-q+t} |\tilde{\mathbf{m}}|_{t,r,T^r} \quad \forall T \in \mathcal{S}_{1,h} \quad \forall \mathbf{m} \in [W^{t,r}(\tilde{\Omega})]^2. \qquad (45) $$

The proof of (45) for $q = 0$, $s = 2$, $t = k + 1$ and $r = 2$ will be given in Section 6 in the framework of an estimate for the interpolation error in the norm $[\cdot]_h$. However, since only Sobolev norms are involved in (45), it is easy to establish that it holds true, thanks to standard estimates (cf. [19, 32]).
(A) The convex case

We pursue the study of the uniform stability of $d_h$ by first assuming that the problem-definition domain is convex. This is because the analysis, though still complex, is simpler in this case.

First of all we introduce a suitable counterpart of space $W_h$, namely,

$$W := \{ \mathbf{q} | \mathbf{q} \in \mathbf{Q} \cap [H^{1/2+\varepsilon}(\Omega)]^2 \},$$

where $\varepsilon$ is a small strictly positive real number. $W$ is equipped with the norm $\|\cdot\|_W := \|\cdot\|_{1/2+\varepsilon} + \|\nabla \cdot\|_0$ (for details on Sobolev spaces of fractional order we refer to [17]).

We next prove

**Proposition 5.4.** For every $v \in V_h$ define a function $\tilde{v}$ in $\Omega$ associated with $v$, as the extension of $v$ by zero in $\bar{\Omega} \setminus \bar{\Omega}_h$. In doing so, it holds

$$\exists \delta_W > 0 \text{ such that } \forall v \in V_h \sup_{\mathbf{p} \in W \setminus \mathbf{0}_W} \frac{(\nabla \cdot \mathbf{p}, \tilde{v})}{\|\mathbf{p}\|_W} \ge \delta_W \|\tilde{v}\|_0. \quad (46)$$

*Proof.* Noting that $\tilde{v} \in L_0^2(\Omega)$ if $length(\Gamma_0) = 0$, we consider a function $\tilde{z}$ defined in $\Omega$ to be the unique solution of the following Poisson problem:

$$\left\{
\begin{array}{ll}
-\Delta \tilde{z} = \tilde{v} & \text{in } \Omega \\
\tilde{z} = 0 & \text{on } \Gamma_0 \text{ if } length(\Gamma_0) > 0 \text{ and } \int_{\Omega} \tilde{z} = 0 \text{ if } length(\Gamma_0) = 0 \\
\partial \tilde{z} / \partial n = 0 & \text{on } \Gamma_1.
\end{array}
\right.
\quad (47)$$

Setting $\tilde{\mathbf{p}} := -\nabla \tilde{z}$ we have $\tilde{\mathbf{p}} \in \mathbf{Q}$ and, according to [22], $\tilde{\mathbf{p}} \in [H^{1/2+\varepsilon}]^2$. Furthermore, there exists a constant $\tilde{C}$ depending only on $\Omega$ and $\Gamma_1$ such that

$$\|\tilde{\mathbf{p}}\|_{1/2+\varepsilon} \le \tilde{C} \|\tilde{v}\|_0. \quad (48)$$

Since $(\nabla \cdot \tilde{\mathbf{p}}, \tilde{v}) = (-\Delta \tilde{z}, \tilde{v}) = \|\tilde{v}\|_0^2$, using (48), (46) is trivially seen to hold with $\delta_W = (1 + \tilde{C})^{-1}$. $\square$

As much as $\Pi_h$ was used in the study of $b_h$, the uniform inf-sup condition for $d_h$ relies on the properties of the interpolation operator $\tilde{\Pi}_h$ defined in (44). Here we take $t = 1/2 + \varepsilon$ and $r = 2$, to have $\tilde{\Pi}_{h}: W \to P_h$ defined as in (44) by replacing $\tilde{\mathbf{m}}$ with $\mathbf{m} \in W$. This is because we can do without extensions in $\bar{\Omega}$, since $\Omega$ is convex and thus the estimate (45) also holds for any $\mathbf{m} \in W$ instead of $\tilde{\mathbf{m}}$.

Let us prove two preliminary results for the interpolation operator $\tilde{\Pi}_h$, namely,

**Proposition 5.5.** There exists a constant $C_{\tilde{\Pi}}$ independent of $h$ such that

$$\left[ \left\| \tilde{\Pi}_h \mathbf{m} \right\| \right]_h \le C_{\tilde{\Pi}} \| \mathbf{m} \|_W \quad \forall \mathbf{m} \in W. \quad (49)$$

*Proof.* First of all, from the same arguments as those applying to the operator $\Pi_h$ (see e.g., [10]), it is easy to figure out that, for a suitable constant $C'_{\tilde{\Pi}}$ independent of $h$, it holds

$$\left\| \tilde{\Pi}_h \mathbf{m} \right\|_{0,h} \le C'_{\tilde{\Pi}} \| \mathbf{m} \|_{0,s,h} \quad \forall \mathbf{m} \in W. \quad (50)$$

Indeed, although zero values are prescribed to $\tilde{\Pi}_h \mathbf{m} \cdot \mathbf{n}$ at some points of $\Gamma_1$, we are not assuming that $\tilde{\Pi}_h \mathbf{m}$ takes up values of $\mathbf{m}$ anywhere in a pointwise sense.
On the other hand, the Rellich–Kondratchev Theorem ensures that $\mathbf{W}$ is continuously embedded in $[L^s(\Omega)]^2$ for $2 < s \le 4$. Hence, from (50) and the convexity of $\Omega$, we immediately infer the existence of another mesh-independent constant $C_{\tilde{\Pi}}^0$ such that

$$ \|\tilde{\Pi}_h \mathbf{m}\|_{0,h} \le C_{\tilde{\Pi}}^0 \|\mathbf{m}\|_{\mathbf{W}} \quad \forall \mathbf{m} \in \mathbf{W}. \qquad (51) $$

Now we turn our attention to $\nabla \cdot \tilde{\Pi}_h \mathbf{m}$, which we denote by $g_h$ for simplicity.

First of all, owing to the definition of $\tilde{\Pi}_h$, $\forall T \in \mathcal{T}_h \setminus S_{1,h}$ we successively have

$$ \left\{ \begin{aligned} \|g_h\|_{0,T}^2 &= \int_T (\nabla \cdot \tilde{\Pi}_h \mathbf{m} - \nabla \cdot \mathbf{m}) \nabla \cdot \tilde{\Pi}_h \mathbf{m} + \int_T \nabla \cdot \tilde{\Pi}_h \mathbf{m} \nabla \cdot \mathbf{m} \\ &= \int_{\partial T} (\tilde{\Pi}_h \mathbf{m} - \mathbf{m}) \cdot \mathbf{n}_T \nabla \cdot \tilde{\Pi}_h \mathbf{m} - \int_T (\tilde{\Pi}_h \mathbf{m} - \mathbf{m}) \cdot \nabla \nabla \cdot \tilde{\Pi}_h \mathbf{m} + \int_T \nabla \cdot \tilde{\Pi}_h \mathbf{m} \nabla \cdot \mathbf{m} \end{aligned} \right. \qquad (52) \\ &= \int_T \nabla \cdot \tilde{\Pi}_h \mathbf{m} \nabla \cdot \mathbf{m}, $$

which yields

$$ \|g_h\|_{0,T} \le \|\nabla \cdot \mathbf{m}\|_{0,T} \quad \forall T \in \mathcal{T}_h \setminus S_{1,h}. \qquad (53) $$

Now taking $T \in S_{1,h}$, again by the definition of $\tilde{\Pi}_h$ we successively have

$$ \|g_h\|_{0,T}^2 = \oint_{\partial T} \tilde{\Pi}_h \mathbf{m} \cdot \mathbf{n}_T g_h - \int_T \tilde{\Pi}_h \mathbf{m} \cdot \nabla g_h, \qquad (54) $$

$$ \|g_h\|_{0,T}^2 = \oint_{\partial T} \mathbf{m} \cdot \mathbf{n}_T g_h - \int_T \mathbf{m} \cdot \nabla g_h + \int_{e_T} [\tilde{\Pi}_h \mathbf{m} - \mathbf{m}] \cdot \mathbf{n}_T g_h, \qquad (55) $$

$$ \|g_h\|_{0,T}^2 = \int_T \nabla \cdot \mathbf{m} g_h - \int_{e_T} \mathbf{m} \cdot \mathbf{n}_T g_h + \int_{e_T} \tilde{\Pi}_h \mathbf{m} \cdot \mathbf{n}_T g_h. \qquad (56) $$

Since $\mathbf{m} \cdot \mathbf{n}$ vanishes on $\Gamma_1$, referring to Figure 1 and using the set $T^\Gamma = T \cup \Delta_T$ we can write

$$ \left\{ \begin{aligned} \|g_h\|_{0,h}^2 &= \sum_{T \in S_{1,h}} [J_{1,T}(\mathbf{m}, g_h) + J_{2,T}(\mathbf{m}, g_h) + J_{3,T}(\mathbf{m}, g_h)] \quad \text{where} \\ J_{1,T}(\mathbf{m}, g_h) &:= \int_{T^\Gamma} \nabla \cdot \mathbf{m} g_h, \\ J_{2,T}(\mathbf{m}, g_h) &:= \int_{\Delta_T} \mathbf{m} \cdot \nabla g_h, \\ J_{3,T}(\mathbf{m}, g_h) &:= \int_{e_T} \tilde{\Pi}_h \mathbf{m} \cdot \mathbf{n}_T g_h. \end{aligned} \right. \qquad (57) $$

By the Cauchy–Schwarz inequality, followed by (11) and Young’s inequality, we obtain for $C_{J1} := C_0^+/C_0^-$

$$ |J_{1,T}(\mathbf{m}, g_h)| \le \frac{\|g_h\|_{0,T}^2}{4} + C_{J1}^2 \|\nabla \cdot \mathbf{m}\|_{0,T^\Gamma}^2. \qquad (58) $$

Now, since $\mathbf{m} \in [L^4(\Omega)]^2$, using Hölder’s inequality with $s = 4$ and $r = 4/3$ we come up with

$$ |J_{2,T}(\mathbf{m}, g_h)| \le \left[ \int_{\Delta_T} |\mathbf{m}|^s \right]^{1/s} \left[ \int_{\Delta_T} |\nabla g_h|^r \right]^{1/r} \le \left[ \int_{\Delta_T} |\mathbf{m}|^s \right]^{1/s} C_0^{1/r} h_T^{3/r} \|\nabla g_h\|_{0, \infty, T^\Gamma}. \qquad (59) $$

Then resorting again to (11) we further have

$$ |J_{2,T}(\mathbf{m}, g_h)| \le C_0^{1/r} C_0^+ \left[ \int_{T^\Gamma} |\mathbf{m}|^s \right]^{1/s} h_T^{3/r-1} \|\nabla g_h\|_{0,T}. \qquad (60) $$

From (60) together with well known inverse inequalities (see e.g., [19]), we infer the existence of a mesh-independent constant $C'_{J2}$ such that

$$ |J_{2,T}(\mathbf{m}, g_h)| \le C'_{J2} \left[ \int_{T^\Gamma} |\mathbf{m}|^s \right]^{1/s} h_T^{3/r-2} \|g_h\|_{0,T}. \qquad (61) $$
As for $J_{3,T}(\mathbf{m}, g_h)$, denoting by $\omega_i$ the weights of the Gauss quadrature formula on $e_T$ we have

$$ 
\begin{aligned}
|J_{3,T}(\mathbf{m}, g_h)| & \le l_T \sum_{i=1}^{k+1} \omega_i |\tilde{\Pi}_h \mathbf{m}|(M_i) \cdot \mathbf{n}_T |g_h(M_i)| \\
& \le h_T \|g_h\|_{0,\infty,T} \sum_{i=1}^{k+1} |\tilde{\Pi}_h \mathbf{m}|(M_i) \cdot \mathbf{n}_T.
\end{aligned} 
\quad (62) $$

Then, recalling the definition of the points $N_i$ on $\Gamma_1$ and the notation $\mathbf{n}_i = \mathbf{n}(N_i)$, (62) yields

$$ |J_{3,T}(\mathbf{m}, g_h)| \le h_T \|g_h\|_{0,\infty,T} \sum_{i=1}^{k+1} |[\tilde{\Pi}_h \mathbf{m}](M_i) \cdot \mathbf{n}_T - [\tilde{\Pi}_h \mathbf{m}](N_i) \cdot \mathbf{n}_i|. \quad (63) $$

Adding and subtracting $[\tilde{\Pi}_h \mathbf{m}](N_i) \cdot \mathbf{n}_T$ inside the absolute value on the right hand side of (63), owing to Propositions 3.1 and 3.2 we easily obtain, for a suitable mesh-independent constant $C'_{J3}$

$$ |J_{3,T}(\mathbf{m}, g_h)| \le C'_{J3} h_T \|g_h\|_{0,\infty,T} \left[ h_T \|\tilde{\Pi}_h \mathbf{m}\|_{0,\infty,T^\Gamma} + h_T^2 \|\nabla \tilde{\Pi}_h \mathbf{m}\|_{0,\infty,T^\Gamma} \right]. \quad (64) $$

Finally, applying standard inverse inequalities (see e.g., [19]) and (11) to (64), followed by Young's inequality, we easily come up with a constant $C_{J3}$ independent of $h$ such that

$$ |J_{3,T}(\mathbf{m}, g_h)| \le \frac{\|g_h\|_{0,T}^2}{4} + C_{J3}^2 \|\tilde{\Pi}_h \mathbf{m}\|_{0,T}^2. \quad (65) $$

Now all that is left to do is to sum up over $S_{1,h}$ the three upper bounds (58), (61) and (65). While on the one hand the first and the third ones are trivially handled, on the other hand we still have some work to do, as far as (61) is concerned.

Using Hölder's inequality with $s = 4$ and $r = 4/3$, together with a trivial upper bound, we obtain

$$ \sum_{T \in S_{1,h}} |J_{2,T}(\mathbf{m}, g_h)| \le C'_{J2} \|\mathbf{m}\|_{0,s} \left[ \sum_{T \in S_{1,h}} h_T^{3-2r} \|g_h\|_{0,T}^r \right]^{1/r} . \quad (66) $$

Next we use again Hölder's inequality in the summation on the right hand side with conjugate exponents $s'$ and $r'$ for $r' = 2/r = 3/2$, to obtain

$$ 
\begin{aligned}
\sum_{T \in S_{1,h}} |J_{2,T}(\mathbf{m}, g_h)| & \le C'_{J2} \|\mathbf{m}\|_{0,s} \left[ \sum_{T \in S_{1,h}} h_T^{(3-2r)s'} \right]^{1/(rs')} \|g_h\|_{0,h} \\
& = C'_{J2} \|\mathbf{m}\|_{0,s} \left[ \sum_{T \in S_{1,h}} h_T \right]^{1/4} \|g_h\|_{0,h}.
\end{aligned} 
\quad (67) $$

Obviously enough, there exists a constant $C_{\Gamma_1}$ independent of $h$ such that $\sum_{T \in S_{1,h}} h_T \le C_{\Gamma_1} \text{length}(\Gamma_1)$. Thus using Young's inequality, for a suitable mesh-independent constant $C_{J2}$, (67) easily yields

$$ \sum_{T \in S_{1,h}} |J_{2,T}(\mathbf{m}, g_h)| \le \frac{\|g_h\|_{0,h}^2}{4} + C_{J2}^2 \|\mathbf{m}\|_{0,4}^2. \quad (68) $$
Summing also up all the terms in (58) and (65) and plugging into (57) the resulting inequalities together with (68), we readily obtain

$$ \sum_{T \in S_{1,h}} \left\| \nabla \cdot \tilde{\Pi}_h \mathbf{m} \right\|_{0,T}^2 \le 2 \{ \max[C_{J1}, C_{J2}, C_{J3}] \}^2 \left\{ \|\mathbf{m}\|_0^2 + \|\mathbf{m}\|_{0,4}^2 + \|\tilde{\Pi}_h \mathbf{m}\|_{0,h}^2 \right\}. \quad (69) $$

Now, taking into account (53), the following upper bound for $\|\nabla \cdot \tilde{\Pi}_h \mathbf{m}\|_{0,h}$ derives from (69):

$$ \left\| \nabla \cdot \tilde{\Pi}_h \mathbf{m} \right\|_{0,h}^2 \le 2 \{ \max[C_{J1}, C_{J2}, C_{J3}] \}^2 \left\{ \|\mathbf{m}\|_0^2 + \|\mathbf{m}\|_{0,4}^2 + \|\tilde{\Pi}_h \mathbf{m}\|_{0,h}^2 \right\} + \|\nabla \cdot \mathbf{m}\|_{0,h}^2. \quad (70) $$

Finally, recalling (51) and the fact that $H^{1/2+\varepsilon}(\Omega)$ is continuously embedded into $L^4(\Omega)$, we come up with a constant $C_{\tilde{\Pi}}^D$ such that

$$ \left\| \nabla \cdot \tilde{\Pi}_h \mathbf{m} \right\|_{0,h} \le C_{\tilde{\Pi}}^D \|\mathbf{m}\|_{\mathbf{W}} \quad \forall \mathbf{m} \in \mathbf{W}, \quad (71) $$

which completes the proof with $C_{\tilde{\Pi}} \sqrt{(C_{\tilde{\Pi}}^0)^2 + (C_{\tilde{\Pi}}^D)^2}$. $\square$

**Proposition 5.6.** There exists a constant $C_\delta$ independent of $h$ such that

$$ d_h(\tilde{\Pi}_h \mathbf{m} - \mathbf{m}, v)_h \le C_\delta h^{1/4} \|\mathbf{m}\|_{\mathbf{W}} \|v\|_{0,h} \quad \forall \mathbf{m} \in \mathbf{W} \text{ and } \forall v \in V_h. \quad (72) $$

*Proof.* From the definition of $\tilde{\Pi}_h$ we may write $\forall v \in V_h$

$$ \begin{aligned} d_h(\tilde{\Pi}_h \mathbf{m} - \mathbf{m}, v) &= \sum_{T \in \mathcal{T}_h} \left\{ \oint_{\partial T} \left[ \tilde{\Pi}_h \mathbf{m} - \mathbf{m} \right] \cdot \mathbf{n}_T v - \left( \tilde{\Pi}_h \mathbf{m} - \mathbf{m}, \nabla v \right)_{0,T} \right\} \\ &= \sum_{T \in \mathcal{S}_{1,h}} \int_{e_T} \left[ \tilde{\Pi}_h \mathbf{m} - \mathbf{m} \right] \cdot \mathbf{n}_T v. \end{aligned} \quad (73) $$

Since $\mathbf{m} \cdot \mathbf{n}$ vanishes identically on $\Gamma_1$, extending $v$ to $\Delta_T$ in the natural way $\forall T \in S_{1,h}$ the following expression derives from (73):

$$ d_h(\tilde{\Pi}_h \mathbf{m} - \mathbf{m}, v) = \sum_{T \in S_{1,h}} \left\{ \int_{e_T} [\tilde{\Pi}_h \mathbf{m} - \mathbf{m}] \cdot \mathbf{n}_T v - \int_{\Gamma_T} [\tilde{\Pi}_h \mathbf{m} - \mathbf{m}] \cdot \mathbf{n} v + \int_{\Gamma_T} \tilde{\Pi}_h \mathbf{m} \cdot \mathbf{n} v \right\}. \quad (74) $$

Now applying the Divergence Theorem in $\Delta_T$, from (74) we obtain rather easily

$$ \left\{ \begin{array}{l} d_h(\tilde{\Pi}_h \mathbf{m} - \mathbf{m}, v) = \sum_{T \in S_{1,h}} [K_{1,T}(\mathbf{m}, v) + K_{2,T}(\mathbf{m}, v) + K_{3,T}(\mathbf{m}, v)] \\ \text{with} \\ K_{1,T}(\mathbf{m}, v) := \int_{\Gamma_T} \tilde{\Pi}_h \mathbf{m} \cdot \mathbf{n} v \\ K_{2,T}(\mathbf{m}, v) := - \int_{\Delta_T} \nabla \cdot (\tilde{\Pi}_h \mathbf{m} - \mathbf{m}) v \\ K_{3,T}(\mathbf{m}, v) := - \int_{\Delta_T} (\tilde{\Pi}_h \mathbf{m} - \mathbf{m}) \cdot \nabla v. \end{array} \right. \quad (75) $$

First of all the functional $K_{1,T}$ can be dealt with by means of Lemma 4.3, that is,

$$ |K_{1,T}(\mathbf{m}, v)| \le C_{\Gamma_1} h_T^{k+1} \left\| \tilde{\Pi}_h \mathbf{m} \right\|_{k+1,T} \|v\|_{0,T}. \quad (76) $$
Now we repeatedly resort to fractional inverse inequalities (cf. [20]) to infer the existence of a constant $C_{1/2}$ independent of $T$ such that

$$|K_{1,T}(\mathbf{m}, v)| \leq C_{1/2} h_T^{1/2} \left\| \tilde{\Pi}_h \mathbf{m} \right\|_{1/2+\varepsilon, T} \|v\|_{0,T}. \quad (77)$$

Owing to the fact that, by construction, to each triangle $T \in S_{1,h}$ exactly $(k+1)(k+2)$ linearly independent DOFs not associated with $\Gamma_T$ are assigned to $\tilde{\Pi}_h$, this operator is invariant as applied to fields in $[P_k(T)]^2$ independently of boundary DOFs. Therefore, on the basis of the well established interpolation theory in $\mathbb{R}^n$, including the case of vector fields (see e.g., [8, 10]), and recalling also approximation results in fractional norms [32], we can assert that there exists another constant $C_\varepsilon$ independent of $T$ such that

$$\left\| \tilde{\Pi}_h \mathbf{m} \right\|_{1/2+\varepsilon, T} \leq C_\varepsilon \| \mathbf{m} \|_{1/2+\varepsilon, T}. \quad (78)$$

Setting $C_{K,1} := C_{1/2}C_\varepsilon$, this immediately yields the following estimate for $K_{1,T}(\mathbf{m}, v)$:

$$|K_{1,T}(\mathbf{m}, v)| \leq C_{K,1} h_T^{1/2} \|\mathbf{m}\|_{1/2+\varepsilon, T} \|v\|_{0,T}. \quad (79)$$

Next we observe that both $\nabla \cdot \tilde{\Pi}_h \mathbf{m}$ and $v$ are polynomials in $P_k(T^\Gamma)$. In doing so we resort to a well known argument (see e.g., [34]), according to which

$$\exists C_\Delta \text{ independent of } T \text{ such that } \|p\|_{0,\Delta_T} \le C_\Delta h_T^{1/2} \|p\|_{0,T} \quad \forall p \in P_k(T^\Gamma) \text{ and } \forall T \in S_{1,h}. \quad (80)$$

Applying (80) to $K_{2,T}(\mathbf{m}, v)$ we readily come up with

$$|K_{2,T}(\mathbf{m}, v)| \leq C_\Delta h_T^{1/2} \left( C_\Delta h_T^{1/2} \left\| \nabla \cdot \tilde{\Pi}_h \mathbf{m} \right\|_{0,T} + \| \nabla \cdot \mathbf{m} \|_{0,T^\Gamma} \right) \|v\|_{0,T}. \quad (81)$$

As for $K_{3,T}(\mathbf{m}, v)$, we first employ Hölder's inequality to obtain

$$|K_{3,T}(\mathbf{m}, v)| \le \left\| \tilde{\Pi}_h \mathbf{m} - \mathbf{m} \right\|_{0,4,T^\Gamma} \|\nabla v\|_{0,4/3,\Delta_T} \quad (82)$$

Next we note that, according to the Rellich–Kondratchev Theorem (cf. [1]) $H^{1/2+\varepsilon}(D)$ is continuously embedded in $L^4(D)$ for all domains $D$ in the aforementioned class. Thus, (45) yields

$$\left\| \tilde{\Pi}_h \mathbf{m} - \mathbf{m} \right\|_{0,4,T^\Gamma} \leq C_{0,4} h_T^\varepsilon |\mathbf{m}|_{1/2+\varepsilon, T^\Gamma}. \quad (83)$$

Taking into account (11), plugging (83) into (82) we come up with

$$|K_{3,T}(\mathbf{m}, v)| \leq C_{0,4} C_0^{3/4} h_T^{9/4+\varepsilon} |\mathbf{m}|_{1/2+\varepsilon, T^\Gamma} \|\nabla v\|_{0,\infty, T^\Gamma} \leq C_0^+ C_{0,4} C_0^{3/4} h_T^{5/4+\varepsilon} |\mathbf{m}|_{1/2+\varepsilon, T^\Gamma} \|\nabla v\|_{0,T}. \quad (84)$$

Owing to (84) and a classical inverse inequality, we guarantee the existence of a constant $C_{1/4}$ independent of $T$ such that

$$|K_{3,T}(\mathbf{m}, v)| \leq C_{1/4} h_T^{1/4} |\mathbf{m}|_{1/2+\varepsilon, T^\Gamma} \|v\|_{0,T}. \quad (85)$$

Now plugging (79), (81) and (85) into (75), we easily obtain for a suitable constant $C'_\Delta$

$$\left| \left( \nabla \cdot \left[ \tilde{\Pi}_h \mathbf{m} - \mathbf{m} \right] \right)_h \right| \leq C'_\Delta h_T^{1/4} \left( \left\| \nabla \cdot \tilde{\Pi}_h \mathbf{m} \right\|_{0,h} + \|\mathbf{m}\|_{\mathbf{W}} + |\mathbf{m}|_{1/2+\varepsilon} \right) \|v\|_{0,h}, \quad (86)$$

Finally, taking into account (71), the result follows. □
As a direct consequence of Propositions 5.5 and 5.6 we have,

**Theorem 5.7.** Provided $h$ is sufficiently small and $\Gamma$ is of the piecewise $C^{2k+3}$-class, there exists a constant $\delta > 0$ independent of $h$, such that,

$$ \forall v \in V_h \quad \sup_{\mathbf{p} \in \mathbf{P}_h \text{ s.t. } [\|\mathbf{p}\|_h = 1]} d_h(\mathbf{p}, v) \geq \delta \|v\|_{0,h}. \quad (87) $$

*Proof.* According to Proposition 5.4, for a given $v \in V_h \ni \tilde{\mathbf{p}} \in \mathbf{W}$ such that

$$ d_h(\tilde{\mathbf{p}}, v) \geq \delta \mathbf{w} \|\tilde{\mathbf{p}}\|_{\mathbf{w}} \|v\|_{0,h}. \quad (88) $$

Combining (72) and (88) we can write

$$ d_h(\tilde{\Pi}_h \tilde{\mathbf{p}}, v) \geq \left( \delta \mathbf{w} - C_\delta h^{1/4} \right) \|\tilde{\mathbf{p}}\|_{\mathbf{w}} \|v\|_{0,h}. \quad (89) $$

Assume that $h^{1/4} < \frac{\delta \mathbf{w}}{2C_\delta}$. In this case (49) allows us to write

$$ \frac{d_h(\tilde{\Pi}_h \tilde{\mathbf{p}}, v)}{[\|\tilde{\Pi}_h \tilde{\mathbf{p}}\|]_h} \geq \frac{\delta \mathbf{w} - C_\delta h^{1/4}}{C_{\tilde{\Pi}}} \|v\|_{0,h}. \quad (90) $$

Thus defining $\mathbf{p}_1 := \frac{\tilde{\Pi}_h \tilde{\mathbf{p}}}{[\|\tilde{\Pi}_h \tilde{\mathbf{p}}\|]_h} \in \mathbf{P}_h$, it holds

$$ d_h(\mathbf{p}_1, v) \geq \frac{\delta \mathbf{w}}{2C_{\tilde{\Pi}}} \|v\|_{0,h}, \quad (91) $$

with $[\|\mathbf{p}_1\|]_h = 1$. Then the result follows with $\delta = \delta_{\mathbf{w}}/(2C_{\tilde{\Pi}})$. $\square$

Finally, taking into account (6), we can prove the main result in this section, namely,

**Theorem 5.8.** Provided $h$ is sufficiently small and $\Omega$ is a convex domain smooth enough, (7) has a unique solution. Moreover there exists a constant $\gamma > 0$ independent of $h$ such that,

$$ \forall (\mathbf{p}; u) \in \mathbf{P}_h \times V_h \quad \sup_{(\mathbf{q}; v) \in \mathbf{Q}_h \times V_h \text{ s.t. } |||(\mathbf{q}; v)|||_h = 1} c_h((\mathbf{p}; u), (\mathbf{q}; v)) \geq \gamma |||(\mathbf{p}; u)|||_h. \quad (92) $$

*Proof.* According to [30] and owing to Theorems 5.1, 5.3 and 5.7, the bilinear form $c_h$ defined by (4) and (5) is weakly coercive over the pair of spaces $\mathbf{P}_h \times U_h$ and $\mathbf{Q}_h \times V_h$, with $U_h = V_h$. This allows us to assert that problem (7) has a unique solution (cf. [16]).

Furthermore, again by [30] (cf. Thm. 3.1), the weak coercivity constant $\gamma$ for bilinear form $c_h$ given by (4) and (5) is also independent of $h$, since the three constants $\alpha, \beta$ and $\delta$ are independent of $h$. Actually, according to [30], $\gamma$ may be expressed by

$$ \left\{ \begin{array}{l} \gamma := \beta\nu/\sqrt{2[(\nu+1)^2 + \beta^2]} \\ \text{with} \\ \nu = \alpha\delta/\sqrt{(1+\delta)^2 + \alpha^2}. \end{array} \right. \quad (93) $$

$\square$
(B) The non convex case

At the price of more or less lengthier technicalities, we next extend the study of the uniform stability of problem (7) for a convex $\Omega$ to non convex domains.

The non convex counterpart of condition (87) will be established by exploiting the so-called extension technique. With this aim, we first observe that it does make sense to define an interpolate in $P_h$ of a field $\tilde{m}$ in a certain analog $\tilde{W}$ of $W$ with respect to $\tilde{\Omega}$, namely, $\tilde{W} := \{\tilde{m}| \tilde{m} \in H(div, \tilde{\Omega}), \tilde{m}|_D \in [H^{1/2+\varepsilon}(D)]^2$, for $D = \Omega$ and $D = \Delta_{\Omega}$ and $\tilde{m} \cdot n|_{\Gamma_1} = 0$ on $\Gamma_1\}$.

We equip $\tilde{W}$ with the norm $\|\cdot\|_{\tilde{W}}$ defined as follows $\forall \tilde{m} \in \tilde{W}$:

$$ \|\tilde{m}\|_{\tilde{W}} := \|\tilde{m}\|_{0,\tilde{\Omega}} + \|\nabla \cdot \tilde{m}\|_{0,\tilde{\Omega}} + \|\tilde{m}|_{\Omega}\|_{1/2+\varepsilon} + \|\tilde{m}|_{\Delta_{\Omega}}\|_{1/2+\varepsilon,\Delta_{\Omega}}. \quad (94) $$

We assume that $\forall m \in W$ there exists a suitable extension $\tilde{m} \in \tilde{W}$ of $m$ to the whole $\tilde{\Omega}$ satisfying for a suitable constant $\tilde{C}_W$

$$ \|\tilde{m}\|_{\tilde{W}} \le \tilde{C}_W \|m\|_{W} \quad \forall m \in W. \quad (95) $$

Once the extension $\tilde{m} \in \tilde{W}$ of $m$ is fixed, the interpolate of $\tilde{m}$ in $P_h$ denoted by $\Pi'_h \tilde{m}$ is the field satisfying

$$ \left\{ \begin{array}{ll} \forall T \in T_h \setminus S_{1,h}, \quad \int_e [\Pi'_h \tilde{m} - \tilde{m}] \cdot n_T v = 0 & \forall v \in P_k(e) \quad \forall \text{ edge } e \subset \partial T; \\ \forall T \in S_{1,h}, \quad \int_e [\Pi'_h \tilde{m} - m] \cdot n_T v = 0 & \forall v \in P_k(e) \quad \forall \text{ edge } e \subset \partial T \setminus e_T. \\ \text{For } k > 0, & \\ \forall T \in T_h \setminus S'_{1,h}, \quad \int_T [\Pi'_h \tilde{m} - \tilde{m}] \cdot o = 0 & \forall o \in [P_{k-1}(T)]^2; \\ \forall T \in S'_{1,h}, \quad \int_{T^{\Gamma'}} [\Pi'_h \tilde{m} - m] \cdot o = 0 & \forall o \in [P_{k-1}(T^{\Gamma'})]^2. \end{array} \right. \quad (96) $$

The existence and uniqueness of $\Pi'_h \tilde{m}$ relies basically upon the same arguments that apply to the interpolation operator $\Pi_h$ from $W_h$ onto $Q_h$, together with Proposition 4.2. In this respect, it is noticeable that for $T \in S'_{0,h}$ the intersection of $e_T$ with $\Gamma$ is necessarily the set of its end-points, by assumption. Therefore the condition on the normal trace of $\Pi'_h \tilde{m}$ over $e_T$ makes sense, since not only $\tilde{m}$ belongs to $[L^s(\Delta_T)]^2$ for some $s > 2$ (cf. [10]), but also the embedding of $H^{1/2+\varepsilon}(\Delta_T)$ into $L^s(\Delta_T)$ is ensured by the Rellich-Kondratchev Theorem [1].

The counterparts of (49) and (72) for operator $\Pi'_h$ will be established hereafter. Before pursuing however, we note that a condition analogous to (46) in Proposition 5.4 holds here, namely,

**Proposition 5.9.**

$$ \left\{ \begin{array}{l} \exists \delta'_{\mathbf{W}} > 0 \text{ such that } \forall T_h \in \mathcal{F}' \text{ and} \\ \forall v \in V_h \ \exists \mathbf{p} \in \mathbf{W} \setminus \mathbf{0}_{\mathbf{W}} \text{ fulfilling } (\nabla \cdot \mathbf{p}, \tilde{v}) \ge \delta'_{\mathbf{W}} \|\mathbf{p}\|_{\mathbf{W}} \|\tilde{v}\|_0, \end{array} \right. \quad (97) $$

where $\tilde{v}$ is the function identical to $v$ in $\Omega_h$ and extended by zero in $\tilde{\Omega} \setminus \tilde{\Omega}_h$.

*Proof.* First we address the case where $\text{length}(\Gamma_0) = 0$.

We observe that, in principle $\tilde{v} \notin L_0^2(\Omega)$. Thus we define a function $v' \in L_0^2(\Omega)$ by $v' = 0$ in $\tilde{\Omega} \setminus \tilde{\Omega}_h$ and $v' := \tilde{v} + c'$ in $\tilde{\Omega}_h$ where $c'$ is a constant. Since $\int_{\Omega_h} \tilde{v} = 0$, we trivially have

$$ 0 = c' \text{area}(\Omega) + \int_{\Omega} \tilde{v} \implies c' = \frac{\int_{\Omega_h \setminus \Omega} \tilde{v}}{\text{area}(\Omega)}. \quad (98) $$

Notice that

$$ \|v'\|_0^2 \ge \|\tilde{v}\|_0^2 - 2|c'| \int_{\Omega} |\tilde{v}| + [c']^2 \text{area}(\Omega). \quad (99) $$

On the other hand we have

$$ |c'| \le \sum_{T \in S'_h} \frac{\int_{\Delta_T} |\tilde{v}|}{\text{area}(\Omega)}. \quad (100) $$
It is clear that

$$ \int_{\Delta_T} |\tilde{v}| \le C_0 h_T^2 l_T \|\tilde{v}\|_{0,\infty,T} \quad \forall T \in S_h'. \qquad (101) $$

Hence, using (11) together with a straightforward argument, we come up with the constant $C_Δ := C_0 C_0^+$ such that

$$ \int_{\Delta_T} |\tilde{v}| \le C_\Delta h_T^{3/2} l_T^{1/2} \|\tilde{v}\|_{0,T^{\Gamma'}} \quad \forall T \in S_h'. \qquad (102) $$

Plugging (102) into (100) we readily obtain

$$ |c'| \le [\text{area}(\Omega)]^{-1} C_\Delta h^{3/2} \sqrt{\sum_{T \in S_h'} l_T} \|\tilde{v}\|_0 \le [\text{area}(\Omega)]^{-1} C_\Delta h^{3/2} [\text{length}(\Gamma)]^{1/2} \|\tilde{v}\|_0, \qquad (103) $$

from which we infer the existence of a constant $C'(\Omega)$ depending only on $\Omega$ such that

$$ |c'| \le C'(\Omega) h^{3/2} \|\tilde{v}\|_0. \qquad (104) $$

Plugging the above upper bound of $|c'|$ into equation (99) we successively obtain

$$ \|v'\|_0^2 \ge \|\tilde{v}\|_0^2 - 2C'(\Omega)h^{3/2}\|\tilde{v}\|_0 \int_{\Omega} |\tilde{v}| + [C'(\Omega)]^2 h^3 \|\tilde{v}\|_0^2 \text{area}(\Omega), \qquad (105) $$

$$ \|v'\|_0^2 \ge \|\tilde{v}\|_0^2 - 2C'(\Omega)h^{3/2}\|\tilde{v}\|_0^2 \sqrt{\text{area}(\Omega)} + [C'(\Omega)]^2 h^3 \|\tilde{v}\|_0^2 \text{area}(\Omega). \qquad (106) $$

Hence setting $C(\Omega) = 2C'(\Omega)\sqrt{\text{area}(\Omega)}$ we conclude that

$$ \|v'\|_0^2 \ge [1 - C(\Omega)h^{3/2}] \|\tilde{v}\|_0^2. \qquad (107) $$

It follows from (107) that there exists a mesh-independent constant $C_{\mathcal{F}'}$ such that $\forall \mathcal{T}_h \in \mathcal{F}'$ it holds

$$ \|v'\|_0 \ge C_{\mathcal{F}'} \|\tilde{v}\|_0 \quad \forall v \in V_h. \qquad (108) $$

To complete this introductory part of the proof, we set $v' = \tilde{v}$ in $\tilde{\Omega}_h$ in case $\text{length}(\Gamma_0) > 0$.

Now we consider a function $z'$ defined in $\Omega$ to be the unique solution of the following Poisson problem:

$$ \begin{cases} -\Delta z' = v' & \text{in } \Omega \\ z' = 0 & \text{on } \Gamma_0 \text{ if } \text{length}(\Gamma_0) > 0 \text{ and } \int_{\Omega} z' = 0 \text{ if } \text{length}(\Gamma_0) = 0 \\ \partial z' / \partial n = 0 & \text{on } \Gamma_1. \end{cases} \qquad (109) $$

Setting $\mathbf{p}' := -\nabla z'$ we have $\mathbf{p}' \in \mathbf{Q}$ and, according to [22], $\mathbf{p}' \in [H^{1/2+\varepsilon}(\Omega)]^2$. Furthermore, there exists a constant $C'_{\mathbf{W}}$ depending only on $\Omega$ and $\Gamma_1$ such that

$$ \|\mathbf{p}'\|_{1/2+\varepsilon} \le C'_{\mathbf{W}} \|v'\|_0. \qquad (110) $$

Since $(\nabla \cdot \mathbf{p}', c') = c' \int_{\Gamma} \mathbf{p}' \cdot \mathbf{n} = 0$ if $\text{length}(\Gamma_0) = 0$, using (108) we obtain

$$ (\nabla \cdot \mathbf{p}', \tilde{v}) = (\nabla \cdot \mathbf{p}', v') = (-\Delta z', v') = \|v'\|_0^2 \ge C_{\mathcal{F}'}^2 \|\tilde{v}\|_0^2. \qquad (111) $$

Thus, choosing $\mathbf{p} = \mathbf{p}'$ and taking into account (110), (97) is trivially seen to hold with $\delta'_{\mathbf{W}} = C_{\mathcal{F}'}^2 / (1 + C'_{\mathbf{W}})$. $\square$

Next we prove two propositions for operator $\Pi'_h$ analogous to Propositions 5.5 and 5.6 for $\Pi_h$.
**Proposition 5.10.** There exists a constant $C'_{\Pi}$ independent of $h$ such that

$$ \|\Pi'_{h} \tilde{\mathbf{m}}\|_{h} \le C'_{\Pi} \|\tilde{\mathbf{m}}\|_{\tilde{\mathbf{W}}} \quad \forall \tilde{\mathbf{m}} \in \tilde{\mathbf{W}}. \qquad (112) $$

*Proof.* First of all, in regard of (96), it is not difficult to figure out that $\Pi'_{h} \tilde{\mathbf{m}} \equiv \mathbf{m}$ in $T^{\Gamma'}$ $\forall \mathbf{m} \in [P_k(T^{\Gamma'})]^2$. In addition to this the interpolate $\Pi'_{h} \tilde{\mathbf{m}}$ restricted to $T^{\Gamma'}$ fulfills all the other conditions specified in classical interpolation theory (cf. [8, 32]) to satisfy for a constant $C'_{\epsilon}$ independent of $h$

$$ \forall q \in [0, 1/2 + \epsilon] \quad \|\Pi'_{h} \tilde{\mathbf{m}} - \mathbf{m}\|_{q, T^{\Gamma'}} \le C'_{\epsilon} h^{1/2+\epsilon-q} \|\mathbf{m}\|_{1/2+\epsilon, T^{\Gamma'}} \quad \forall T \in \mathcal{T}_h \text{ and } \forall \mathbf{m} \in \mathbf{W}. \qquad (113) $$

It trivially follows from (113) that

$$ \|\Pi'_{h} \tilde{\mathbf{m}}\|'_{0,h} \le (1 + C'_{\epsilon} h^{1/2+\epsilon}) \|\mathbf{m}\|_{1/2+\epsilon} \quad \forall \mathbf{m} \in \mathbf{W}. \qquad (114) $$

On the other hand, recalling (11), we have

$$ \|\Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^{2} \le (C_{0}^{-})^{-2} (C_{0}^{+})^{2} \|\Pi'_{h} \tilde{\mathbf{m}}\|_{0,T^{\Gamma'}}^{2} \quad \forall T \in \mathcal{T}_{h}. \qquad (115) $$

Therefore, summing up both sides of (115) over the mesh, taking into account (114) and assuming that $h < 1$, the following bound is seen to hold with $C'_{\Pi'} = (1 + C'_{\epsilon})C_{0}^{+}/C_{0}^{-}$:

$$ \|\Pi'_{h} \tilde{\mathbf{m}}\|_{0,h} \le C'_{\Pi'} \|\mathbf{m}\|_{\mathbf{W}} \quad \forall \mathbf{m} \in \mathbf{W}. \qquad (116) $$

As for the divergence term, first we note that the case of $\|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}$ for $T \in \mathcal{T}_h \setminus S'_h$ was already addressed in Section 5. Indeed, from (53) together with (69) and recalling (71), we easily infer that

$$ \sum_{T \in \mathcal{T}_h \setminus S'_h} \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^2 \le (C_{\Pi}^D)^2 \|\mathbf{m}\|_{\mathbf{W}}^2 \quad \forall \mathbf{m} \in \mathbf{W}. \qquad (117) $$

Now for triangles in $S'_h$ we first have

$$ \left\{ \begin{array}{ll} \forall T \in S'_{1,h} & \\ \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^2 = J_{1,T}^1(\mathbf{m}) + J_{2,T}^1(\mathbf{m}) + J_{3,T}^1(\mathbf{m}) & \\ \text{with} & \\ J_{1,T}^1(\mathbf{m}) := \int_{T^{\Gamma'}} \nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}(\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}} - \nabla \cdot \mathbf{m}); & (118) \\ J_{2,T}^1(\mathbf{m}) := \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0, \Delta_T}^2; & \\ J_{3,T}^1(\mathbf{m}) := \int_{T^{\Gamma'}} \nabla \cdot \Pi'_{h} \tilde{\mathbf{m}} \nabla \cdot \mathbf{m}. & \end{array} \right. $$

From the definition of the operator $\Pi'_{h}$, integrating by parts we obtain after straightforward calculations

$$ J_{1,T}^{1}(\mathbf{m}) = (\Pi'_{h} \tilde{\mathbf{m}} \cdot \mathbf{n}, \nabla \cdot \Pi'_{h} \tilde{\mathbf{m}})_{0, \Gamma_T}. \qquad (119) $$

Then Lemma 4.3 yields

$$ |J_{1,T}^{1}(\mathbf{m})| \le C'_{\Gamma_1} h_T^{k+1} \|\Pi'_{h} \tilde{\mathbf{m}}\|_{k+1, T^{\Gamma'}} \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0, T}, \qquad (120) $$

from which we obtain in turn for a typical inverse inequality constant $C'_{1/2}$ (cf. [20])

$$ |J_{1,T}^{1}(\mathbf{m})| \le C'_{1/2} h_T^{1/2} \|\Pi'_{h} \tilde{\mathbf{m}}\|_{1/2+\epsilon, T^{\Gamma'}} \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}. \qquad (121) $$

Further using (113) with $q = 1/2 + \epsilon$ and setting $C'_{J,1} = C'_{1/2}(1 + C'_{\epsilon})$, from (121) we infer that

$$ |J_{1,T}^{1}(\mathbf{m})| \le C'_{J,1} h_T^{1/2} \|\mathbf{m}\|_{1/2+\epsilon, T^{\Gamma'}} \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}. \qquad (122) $$
Then using Young's inequality we come up with

$$|J_{1,T}^{1}(\mathbf{m})| \leq (C_{J,1}^{1})^{2} h_{T} \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^{2}/2 + \|\mathbf{m}\|_{1/2+\varepsilon, T^{\tau'}}^{2}/2. \quad (123)$$

As for $J_{2,T}^1(\mathbf{m})$, by a classical inverse inequality, we obtain for a suitable constant $C_{J,2}^1$

$$|J_{2,T}^{1}(\mathbf{m})| \leq C_0 h_T^3 \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,\infty,T}^2 \leq C_{J,2}^1 h_T \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^2. \quad (124)$$

Finally, using Young's inequality, we easily obtain for $J_{3,T}^1(\mathbf{m})$

$$|J_{3,T}^{1}(\mathbf{m})| \leq \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^{2}/2 + \|\nabla \cdot \mathbf{m}\|_{0,T^{\tau'}}^{2}/2 \leq \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^{2}/2 + \|\nabla \cdot \tilde{\mathbf{m}}\|_{0,T}^{2}/2. \quad (125)$$

Collecting (123), (124) and (125), (118) is readily seen to yield

$$\forall T \in S'_{1,h} \quad \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^{2} \leq \left\{1 - 2h_T \left[C_{J,2}^{1} + (C_{J,1}^{1})^{2}\right]\right\}^{-1} \left(\|\mathbf{m}\|_{1/2+\varepsilon, T^{\tau'}}^{2} + \|\nabla \cdot \tilde{\mathbf{m}}\|_{0,T}^{2}\right). \quad (126)$$

Thus, as long as $h$ is sufficiently small, there exists a mesh-independent constant $C_J^1$ such that

$$\forall T \in S'_{1,h} \quad \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^{2} \leq C_{J}^{1} \left( \|\mathbf{m}\|_{1/2+\varepsilon, T^{\tau'}}^{2} + \|\nabla \cdot \tilde{\mathbf{m}}\|_{0,T}^{2} \right). \quad (127)$$

Next we consider triangles in $S'_{0,h}$. We have

$$\left\{
\begin{array}{ll}
\forall T \in S'_{0,h} & \\
\|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^{2} = J_{1,T}^{0}(\mathbf{m}) + J_{2,T}^{0}(\mathbf{m}) & \\
\text{with} & \\
J_{1,T}^{0}(\mathbf{m}) := \int_{T}(\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}} - \nabla \cdot \tilde{\mathbf{m}})\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}; & \\
J_{2,T}^{0}(\mathbf{m}) := \int_{T}\nabla \cdot \tilde{\mathbf{m}}\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}. &
\end{array}
\right.
\quad (128)$$

Taking into account the definition $\Pi'_h$ we can assert that

$$J_{1,T}^{0}(\mathbf{m}) = \int_{\partial T} (\Pi'_{h} \tilde{\mathbf{m}} - \tilde{\mathbf{m}}) \cdot \mathbf{n}_{T} \nabla \cdot \Pi'_{h} \tilde{\mathbf{m}} - \int_{T} (\Pi'_{h} \tilde{\mathbf{m}} - \tilde{\mathbf{m}}) \cdot \nabla \nabla \cdot \Pi'_{h} \tilde{\mathbf{m}} = 0. \quad (129)$$

On the other hand, we have

$$|J_{2,T}^{0}(\mathbf{m})| \leq \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T}^{2}/2 + \|\nabla \cdot \tilde{\mathbf{m}}\|_{0,T}^{2}. \quad (130)$$

Combining (128)–(130) we obtain

$$\forall T \in S'_{0,h} \quad \|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,T} \leq \|\nabla \cdot \tilde{\mathbf{m}}\|_{0,T}. \quad (131)$$

Finally, as an immediate consequence of (117), (127) and (131), it holds for a suitable mesh-independent constant $C_{\Pi'}^D$

$$\|\nabla \cdot \Pi'_{h} \tilde{\mathbf{m}}\|_{0,h} \leq C_{\Pi'}^{D} \|\tilde{\mathbf{m}}\|_{\tilde{\mathbf{W}}} \quad \forall \mathbf{m} \in \mathbf{W}. \quad (132)$$

The addition of (132) and (116) yields (112) with $C'_{\Pi} = \sqrt{(C_{\Pi'}^0)^2 + (C_{\Pi'}^D)^2}$. $\square$

**Proposition 5.11.** There exists a constant $C'_\delta$ independent of $h$ such that $\forall \mathbf{m} \in \mathbf{W}$ and $\forall v \in V_h$

$$d_h(\Pi'_h \tilde{\mathbf{m}}, v) \geq d_h(\tilde{\mathbf{m}}, v) - \sum_{T \in S'_{1,h}} |(\nabla \cdot \tilde{\mathbf{m}}, v)_{\Delta_T}| - C'_\delta h^{1/4} \|\mathbf{m}\|_{\mathbf{W}} \|v\|_{0,h}. \quad (133)$$
Proof. Using the extension $\tilde{\mathbf{m}} \in \tilde{\mathbf{W}}$ of $\mathbf{m}$ satisfying (95), first we note that by construction we have $\forall \mathbf{m} \in \mathbf{W}$ and $\forall v \in V_h$

$$ \sum_{T \in \mathcal{T}_h \setminus S_{1,h}} (\nabla \cdot \Pi'_h \tilde{\mathbf{m}}, v)_T = \sum_{T \in \mathcal{T}_h \setminus S_{1,h}} (\nabla \cdot \tilde{\mathbf{m}}, v)_T. \quad (134) $$

Furthermore, on the basis of Proposition 5.6, we may write $\forall \mathbf{m} \in \mathbf{W}$ and $\forall v \in V_h$

$$ \sum_{T \in S_{1,h} \setminus S'_{1,h}} (\nabla \cdot \Pi'_h \tilde{\mathbf{m}}, v)_T \geq \sum_{T \in S_{1,h} \setminus S'_{1,h}} (\nabla \cdot \tilde{\mathbf{m}}, v)_T - C_\delta h^{1/4} \|\mathbf{m}\|_{\mathbf{w}} \|v\|_{0,h}. \quad (135) $$

As for $T \in S'_{1,h}$, we may write $\forall \mathbf{m} \in \mathbf{W}$ and $\forall v \in V_h$

$$ \left\{ \begin{array}{ll} \forall T \in S'_{1,h} (\nabla \cdot \Pi'_h \tilde{\mathbf{m}}, v)_T = K'_{1,T}(\mathbf{m}, v) + K'_{2,T}(\mathbf{m}, v) + K'_{3,T}(\mathbf{m}, v) & \\ \text{with} & \\ K'_{1,T}(\mathbf{m}, v) = (\nabla \cdot \Pi'_h \tilde{\mathbf{m}}, v)_{\Delta_T} & (136) \\ K'_{2,T}(\mathbf{m}, v) = (\nabla \cdot [\Pi'_h \tilde{\mathbf{m}} - \mathbf{m}], v)_{T^{\Gamma'}} & \\ K'_{3,T}(\mathbf{m}, v) = (\nabla \cdot \mathbf{m}, v)_{T^{\Gamma'}} & \end{array} \right. $$

Using arguments repeatedly exploited in this work we may write for a constant $C'_{K,1}$ independent of $T$

$$ |K'_{1,T}(\mathbf{m}, v)| \leq C'_{K,1} h_T \|\nabla \cdot \Pi'_h \tilde{\mathbf{m}}\|_{0,T} \|v\|_{0,T} \quad \forall T \in S'_{1,h}. \quad (137) $$

Moreover by construction we have $K'_{2,T}(\mathbf{m}, v) = (\Pi'_h \tilde{\mathbf{m}}, v)_{\Gamma_T}$. Thus Lemma 4.3 and (11) yield

$$ |K'_{2,T}(\mathbf{m}, v)| = \left| \int_{\Gamma_T} \Pi'_h \tilde{\mathbf{m}} v \, ds \right| \leq C'_{\Gamma_1} h_T^{k+1} \|\Pi'_h \tilde{\mathbf{m}}\|_{k+1, T^{\Gamma'}} \|v\|_{0, T^{\Gamma'}} \quad \forall T \in S'_{1,h}, \quad (138) $$

where $C'_{\Gamma_1}$ is a constant independent of $T$. Then mimicking (121) and taking $q = 1/2$ in (113), we can assert the existence of another mesh-independent constant $C'_{K,2}$ such that

$$ |K'_{2,T}(\mathbf{m}, v)| \leq C'_{K,2} h_T^{1/2} \|\mathbf{m}\|_{1/2+\varepsilon, T^{\Gamma'}} \|v\|_{0,T} \quad \forall T \in S'_{1,h}. \quad (139) $$

Putting together (137), (139) and (136), $\forall \mathbf{m} \in \mathbf{W}$ and $\forall v \in V_h$ we have

$$ \sum_{T \in S'_{1,h}} (\nabla \cdot \Pi'_h \tilde{\mathbf{m}}, v)_T \geq \sum_{T \in S'_{1,h}} \left[ (\nabla \cdot \mathbf{m}, v)_{T^{\Gamma'}} - \left( C'_{K,1} h_T \|\nabla \cdot \Pi'_h \tilde{\mathbf{m}}\|_{0,T} + C'_{K,2} h_T^{1/2} \|\mathbf{m}\|_{1/2+\varepsilon, T^{\Gamma'}} \right) \|v\|_{0,T} \right]. \quad (140) $$

Thus, taking into account (134), (135) and (140) and setting $C'_K := C'_{K,2} + C_\delta$, since $h < 1$ we obtain $\forall \mathbf{m} \in \mathbf{W}$ and $\forall v \in V_h$

$$ \left\{ \begin{array}{ll} d_h(\Pi'_h \tilde{\mathbf{m}}, v) \geq d_h(\tilde{\mathbf{m}}, v) - \sum_{T \in S'_{1,h}} |(\nabla \cdot \tilde{\mathbf{m}}, v)_{\Delta_T}| & \\ -[C'_K h^{1/4} \|\mathbf{m}\|_{\mathbf{w}} + C'_{K,1} h \|\nabla \cdot \Pi'_h \tilde{\mathbf{m}}\|_{0,h}] \|v\|_{0,h}. & \end{array} \right. \quad (141) $$

Finally, recalling (132) and setting $C'_\delta = C'_K + \tilde{C}_w C'_{\Pi'} C'_{K,1}$, (133) readily derives from (141). $\square$

Propositions 5.9–5.11 are all that is needed to establish the validity of an inf-sup condition analogous to (87) in the non convex case. More precisely we have

**Theorem 5.12.** If $\Gamma$ is of the piecewise $C^{2k+3}$-class, there exists a constant $\delta' > 0$ independent of $h$, such that $\forall T_h \in \mathcal{F}'$ it holds

$$ \forall v \in V_h \quad \sup_{\mathbf{p} \in \mathbf{P}_{k,h}} \quad d_h(\mathbf{p}, v)_h \geq \delta' \|v\|_{0,h}. \quad (142) $$
*Proof.* Recalling Proposition 5.9, given $v \in V_h$, let $\mathbf{p}$ be the field in $\mathbf{W}$ fulfilling (97). The first step to prove (142) is to construct a suitable extension $\tilde{\mathbf{p}} \in \tilde{\mathbf{W}}$ of $\mathbf{p}$ to $\tilde{\Omega}$ satisfying (95). With this aim we pose the following Poisson problem in the outer domain $\Delta_{\Omega}$, in terms of the solution $z'$ of (109):

$$ egin{cases} 
\text{Find } z_{\Delta} \text{ such that} \\ 
-\Delta z_{\Delta} = 0 & \text{in } \Delta_{\Omega} \\ 
z_{\Delta} = 0 & \text{on } \tilde{\Gamma} \\ 
-\partial z_{\Delta}/\partial n = \partial z'/\partial n & \text{on } \Gamma. 
\end{cases} 
\quad (143) $$

Now we set $\tilde{\mathbf{p}} = \mathbf{p}(-\nabla z')$ in $\Omega$ and $\tilde{\mathbf{p}} = -\nabla z_{\Delta}$ in $\Delta_{\Omega}$. It is clear that $\tilde{\mathbf{p}} \in \mathbf{H}(\mathrm{div}, \tilde{\Omega})$, owing to the continuity of the normal trace of $\tilde{\mathbf{p}}$ across $\Gamma$ by construction. Moreover, since $\partial z'/\partial n \in H^{\varepsilon}(\Gamma)$, we have $z_{\Delta} \in H^{3/2+\varepsilon}(\Delta_{\Omega})$, and the harmonic function $z_{\Delta}$ fulfills for a constant $C_{\mathrm{out}}$

$$ \|z_{\Delta}\|_{3/2+\varepsilon, \Delta_{\Omega}} \leq C_{\mathrm{out}} \|\partial z'/\partial n\|_{\varepsilon, \Gamma}. \quad (144) $$

Thus, by the Trace Theorem in $\Omega$, there is another constant $C_{\mathrm{in}}$ such that

$$ \|z_{\Delta}\|_{3/2+\varepsilon, \Delta_{\Omega}} \leq C_{\mathrm{in}} \|z'\|_{3/2+\varepsilon}. \quad (145) $$

Using classical Friedrichs-Poincaré inequalities (cf. [18]), we can assert the existence of a constant $C_{\Delta}$ depending only on $\Omega$ and $\Delta_{\Omega}$ such that

$$ \|\tilde{\mathbf{p}}\|_{1/2+\varepsilon, \Delta_{\Omega}} \leq C_{\Delta} \|\mathbf{p}'\|_{1/2+\varepsilon}. \quad (146) $$

Since $\nabla \cdot \tilde{\mathbf{p}} = 0$ in $\Delta_{\Omega}$, the norm of $\tilde{\mathbf{p}}$ in $\mathbf{H}(\mathrm{div}, \Delta_{\Omega})$ is bounded above by $\|\tilde{\mathbf{p}}\|_{1/2+\varepsilon, \Delta_{\Omega}}$. Hence for $\tilde{C}_{\mathbf{W}} := (1 + C_{\Delta}^2)^{1/2}$ it holds

$$ \|\tilde{\mathbf{p}}\|_{\tilde{\mathbf{W}}} \leq \tilde{C}_{\mathbf{W}} \|\mathbf{p}\|_{\mathbf{W}}. \quad (147) $$

Now for an arbitrary function in $v \in V_h$, $\mathbf{p}$ being the field in $\mathbf{W}$ fulfilling (97) and $\tilde{\mathbf{p}} \in \tilde{\mathbf{W}}$ being the extension of $\mathbf{p}$ fulfilling (146), we have

$$ d_h(\Pi'_h \tilde{\mathbf{p}}, v) = d_h(\Pi'_h \tilde{\mathbf{p}} - \tilde{\mathbf{p}}, v) + \sum_{T \in S'_h} (\nabla \cdot \tilde{\mathbf{p}}, v)_{\Delta_T} + (\nabla \cdot \mathbf{p}, \tilde{v}). \quad (148) $$

Since $\nabla \cdot \tilde{\mathbf{p}} = 0$ in $\Delta_{\Omega}$, and thus in $\Delta_T$ too $\forall T \in S'_h$, taking into account (133), (148) yields

$$ d_h(\Pi'_h \tilde{\mathbf{p}}, v) \geq (\nabla \cdot \mathbf{p}, \tilde{v}) - C'_{\delta} h^{1/4} \|\tilde{\mathbf{p}}\|_{\tilde{\mathbf{W}}} \|v\|_{0,h}. \quad (149) $$

Furthermore, using (97) and (147), from (149) we obtain

$$ d_h(\Pi'_h \tilde{\mathbf{p}}, v) \geq \delta'_{\mathbf{W}} \|\mathbf{p}\|_{\mathbf{W}} \|\tilde{v}\|_0 - C'_{\delta} \tilde{C}_{\mathbf{W}} h^{1/4} \|\mathbf{p}\|_{\mathbf{W}} \|v\|_{0,h}. \quad (150) $$

On the other hand, by a straightforward argument similar to those already exploited in this work, we can assert that, provided $h$ is small enough, there exists a mesh-independent constant $C'_{V}$ such that

$$ \|\tilde{v}\|_{0}^{2} = \|v\|_{0,h}^{2} - \sum_{T \in S'_{h}} \int_{\Delta_{T}} v^{2} \geq \|v\|_{0,h}^{2} - C_{0}h_{T}^{3} \sum_{T \in S'_{h}} \|v\|_{0,\infty,T}^{2} \geq (C'_{V})^{2} \|v\|_{0,h}^{2}. \quad (151) $$

Plugging (151) into (150), it follows that

$$ d_h(\Pi'_h \tilde{\mathbf{p}}, v) \geq (\delta'_{\mathbf{W}} C'_{V} - C'_{\delta} \tilde{C}_{\mathbf{W}} h^{1/4}) \|\mathbf{p}\|_{\mathbf{W}} \|v\|_{0,h}. \quad (152) $$
On the other hand, recalling (112) and (147) we have

$$[|\Pi'_{h}\tilde{\mathbf{p}}|]_{h} \leq C'_{\Pi} \|\tilde{\mathbf{p}}\|_{\bar{\mathbf{w}}} \leq C'_{\Pi} \tilde{C}_{W} \|\mathbf{p}\|_{\mathbf{w}}, \quad (153)$$

and hence from (152) we obtain

$$d_h(\Pi'_h \tilde{\mathbf{p}}, v) \geq \frac{\delta'_W C'_V - C'_\delta \tilde{C}_W h^{1/4}}{C'_\Pi \tilde{C}_W} [|\Pi'_h \tilde{\mathbf{p}}|]_h \|v\|_{0,h}. \quad (154)$$

Thus, provided $h$ is small enough and in particular $h^{1/4} \le \frac{\delta'_W C'_V}{2C'_\delta \tilde{C}_W}$, (142) immediate follows from (154) with $\delta' = \delta'_W C'_V / (2C'_\Pi \tilde{C}_W)$. $\square$

To conclude we have

**Theorem 5.13.** Assume that $\Omega$ is a non convex domain of the piecewise $C^{2k+3}$-class. Provided $h$ is sufficiently small problem (7) has a unique solution. Moreover we can assert that a uniform stability condition of the same type as (92) holds true for this problem.

*Proof.* Recalling that the validity of both (42) and (36) is guaranteed, irrespective of the fact that $\Omega$ be convex or not, from (142) we can assert that the approximate problem (7) has a unique solution. Furthermore recalling (93), a uniform stability condition of the type (92) does hold with a constant $\gamma'$ instead of $\gamma$, given by

$$\left\{ \begin{array}{ll} \gamma' := \beta\nu' / \sqrt{2[(\nu' + 1)^2 + \beta^2]} & \text{with} \\ \nu' = \alpha\delta' / \sqrt{(1 + \delta')^2 + \alpha^2}. & \end{array} \right. \quad (155)$$

$\square$

## 6. INTERPOLATION ERRORS

Generalizing classical estimates for interpolation errors applying to the space $\mathcal{Q}_h$ (cf. [25]), we establish here new ones for the space $\mathbf{P}_h$.

In this section the notations $D^1 w$ and $D^2 w$ will eventually be replaced by $\nabla w$ and $H(w)$ for the gradient and the Hessian of $w$, respectively. Furthermore in the sequel several arguments will rely upon classical inverse inequalities for polynomials defined in $T \in \mathcal{T}_h$ (cf. [19,36]) and their extensions to neighboring sets. Besides (10) and (11) we shall use the following one:

There exists a constant $C_I$ depending only on $k$ and the shape regularity of $\mathcal{T}_h$ (cf. [8], Chap. 4, Sect. 4) such that for any set $T'$ satisfying $T^{\Gamma'} \subseteq T' \subseteq T^{\Gamma}$, it holds for $1 \le j \le k$:

$$\|D^j w\|_{0,T'} \le C_I h_T^{-1} \|D^{j-1} w\|_{0,T'} \quad \forall w \in \mathcal{P}_k(T'). \quad (156)$$

Now let $\varepsilon > 0$ be a small real number. We define

$$\mathcal{M}^k(\Omega) := \left\{ \mathbf{m} \mid \mathbf{m} \in \left[ H^{k+1+\varepsilon(1-\min[k,1])}(\Omega) \right]^2, \nabla \cdot \mathbf{m} \in H^{k+1}(\Omega), \text{ and } [\mathbf{m} \cdot \mathbf{n}]_{|\Gamma_1} = 0 \right\}.$$

We shall also use extensions of fields in $\mathcal{M}^k(\Omega)$ to a domain $\tilde{\Omega}$ containing both $\Omega$ and all the sets $\Omega_h$ for $\mathcal{T}_h \in \mathcal{F}'$. More precisely, here we consider $\tilde{\Omega}$ to be a domain whose boundary $\tilde{\Gamma}$ is such that the minimum of the distance from all of its points to $\Gamma$ is not less than $h_{\max}^2$, where $h_{\max} = 2 \max_{\mathcal{T}_h \in \mathcal{F}'} [h]$. Let

$$\mathcal{M}^k(\tilde{\Omega}) := \left\{ \tilde{\mathbf{m}} \mid \tilde{\mathbf{m}} \in \left[ H^{k+1+\varepsilon(1-\min[k,1])}(\tilde{\Omega}) \right]^2, \nabla \cdot \tilde{\mathbf{m}} \in H^{k+1}(\tilde{\Omega}), \text{ and } [\tilde{\mathbf{m}} \cdot \mathbf{n}]_{|\Gamma_1} = 0 \right\}.$$
With every $\mathbf{m} \in \mathcal{M}^k(\Omega)$ we associate its extension $\tilde{\mathbf{m}}$ to $\tilde{\Omega} \setminus \Omega$ belonging to $\mathcal{M}^k(\tilde{\Omega})$ constructed, for instance, as advocated by Stein et al. [33] (i.e., using polyharmonic functions).

Now, since the restriction to $\Omega$ of any field in $\mathcal{M}^k(\tilde{\Omega})$ necessarily belongs to $\mathbf{W}$ for all $k \ge 0$, given $\tilde{\mathbf{m}}$ in $\mathcal{M}^k(\tilde{\Omega})$ we may define its interpolate in $\mathbf{P}_h$ inside $\tilde{\Omega}_h$ in the same manner as (44) for fields in $\mathbf{W}$ and a convex $\Omega$. In short, we still denote such an interpolation operator from $\mathcal{M}^k(\tilde{\Omega})$ onto $\mathbf{P}_h$ by $\tilde{\Pi}_h$, which is constructed as prescribed in (44).

Notice that here we have $\tilde{\Pi}_h \tilde{\mathbf{m}} \in \mathbf{H}(\mathrm{div}, \tilde{\Omega}_h)$. It is also noteworthy that every $\tilde{\mathbf{m}} \in \mathcal{M}^k(\tilde{\Omega})$ belongs to $[H^2(\tilde{\Omega})]^2$ for every $k \ge 1$ and hence, from the Sobolev Embedding Theorem, it makes sense to prescribe $[\tilde{\Pi}_h \tilde{\mathbf{m}} \cdot \mathbf{n}](P) = [\tilde{\mathbf{m}} \cdot \mathbf{n}](P)$ at any point $P \in \Gamma$ (cf. [1]). Notice that, in case $k = 0$, such a prescription is also admissible, thanks to the assumption that $\mathbf{m} \in [H^{1+\epsilon}(\Omega)]^2$. Moreover, since $\forall \tilde{\mathbf{m}} \in \mathcal{M}^k(\tilde{\Omega})$ and $\forall k \ge 1$ we have $\nabla \cdot \tilde{\mathbf{m}} \in H^2(\Omega)$, for the same reason it is legitimate to associate with $\nabla \cdot \tilde{\mathbf{m}}$ a function $d_h(\tilde{\mathbf{m}}) \in L^2(\Omega_h)$, whose restriction to every $T \in \mathcal{T}_h$ is the standard Lagrange interpolate of $[\nabla \cdot \tilde{\mathbf{m}}]|_T$ in $P_k(T)$ at the $(k+2)(k+1)/2$ lattice points of $T$ (cf. [37]). On the other hand, for $k = 0$, $d_h(\tilde{\mathbf{m}})$ is defined to be the mean value in $T$ of $[\nabla \cdot \tilde{\mathbf{m}}]|_T \forall T \in \mathcal{T}_h$. For every $T \in \mathcal{S}_h$ such that $\Delta_T$ is not contained in $T$, we extend $d_h(\tilde{\mathbf{m}})$ to $\Delta_T$ in the natural manner.

Now, using the affine Piola transformation for $T \in \mathcal{T}_h$ (see e.g., [14]) from $\mathbf{H}(\mathrm{div}, T)$ onto $\mathbf{H}(\mathrm{div}, \hat{T})$ where $\hat{T}$ is the standard unit master triangle, we may extend the estimate given in Proposition 3.6 of Chapter III from [10] to the case of Sobolev spaces of order greater than one. In a more standard way, as in [13] for the second estimate, this leads to the existence of two constants $\tilde{C}_{k,\tilde{\Omega}}$ and $\tilde{C}_{k,d}$ such that, $\forall T \in \mathcal{T}_h \setminus \mathcal{S}_{1,h}$, it holds $\forall \tilde{\mathbf{m}} \in \mathcal{M}^k(\tilde{\Omega})$

$$ \left\{ \begin{array}{ll} \| D^j [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}] \|_{0,T} & \le \tilde{C}_{k,\tilde{\Omega}} h^{k+1-j} |\tilde{\mathbf{m}}|_{k+1,T} & \text{for } j = 0, 1, \dots, k+1 \\ \| \nabla \cdot \tilde{\mathbf{m}} - d_h(\tilde{\mathbf{m}}) \|_{0,T} & \le \tilde{C}_{k,d} h^{k+1} |\nabla \cdot \tilde{\mathbf{m}}|_{k+1,T}. & \end{array} \right. \qquad (157) $$

The estimates (157) are completed by means of the following technical lemma:

**Lemma 6.1.** Provided $h$ is sufficiently small, there exist two mesh-independent constants $\tilde{C}_{k,\tilde{\Omega}}$ and $\tilde{C}_{k,d}$ such that $\forall \tilde{\mathbf{m}} \in \mathcal{M}^k(\tilde{\Omega})$ fulfilling $\tilde{\mathbf{m}} \cdot \mathbf{n}|_{\Gamma_1} = 0$, it holds $\forall T \in \mathcal{S}_{1,h}$

$$ \left\{ \begin{array}{ll} \| D^j [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}] \|_{0,T} & \le \tilde{C}_{k,\tilde{\Omega}} h^{k+1-j} |\tilde{\mathbf{m}}|_{k+1+\varepsilon(1-\min[k,1]),T^\Gamma} & \text{for } j = 0, 1, \dots, k+1 \\ \| \nabla \cdot \tilde{\mathbf{m}} - d_h(\tilde{\mathbf{m}}) \|_{0,T} & \le \tilde{C}_{k,d} h^{k+1} |\nabla \cdot \tilde{\mathbf{m}}|_{k+1,T^\Gamma}. & \end{array} \right. \qquad (158) $$

*Proof.* First of all we point out that the second estimate in (158) is standard, for it results from well known arguments applying to affine transformations between domains with aspect ratio uniformly bounded away from zero, such as $T^\Gamma$ (cf. [8] and the proof Lemma 4.1 of [27]).

As for the first estimate, we use the inverse of the affine Piola transformation $\mathcal{P}_{T^\Gamma}$ from $[L^2(T^\Gamma)]^2$ onto $[L^2(\hat{T}^\Gamma)]^2$, where $\hat{T}^\Gamma$ is a triangle having two perpendicular straight edges with unit length and at most one curved edge. Then the remainder of the proof follows the same principles as in Lemma 4.1 of [27], by extending the approximation results given in Proposition 3.6 of Chapter III of [10] to hilbertian Sobolev spaces of order greater than one and using standard interpolation theory (cf. [13, 32]). $\square$

An important consequence of Lemma 6.1 is

**Theorem 6.2.** Provided $h$ is sufficiently small there exists a mesh-independent constant $C_H$ such that $\forall \tilde{\mathbf{m}} \in \mathcal{M}^k(\tilde{\Omega})$ it holds

$$ \left[ |\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}| \right]_h \le C_H h^{k+1} \left[ \|\tilde{\mathbf{m}}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}^2 + |\nabla \cdot \tilde{\mathbf{m}}|_{k+1,\tilde{\Omega}}^2 \right]^{1/2}. \qquad (159) $$
Proof. First of all, taking $j = 0$ in the first equation of both (157) and (158) and setting $C_{k,\bar{\Omega}} = \max[\tilde{C}_{k,\bar{\Omega}}, \tilde{C}_{k,\bar{\Omega}}]$, we immediately infer that

$$ \|\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}\|_{0,h} \le C_{k,\bar{\Omega}} h^{k+1} |\tilde{\mathbf{m}}|_{k+1+\varepsilon(1-\min[k,1]),\bar{\Omega}}. \quad (160) $$

Next we estimate the term $\|\nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}]\|_{0,h}$.

From standard results for the space $RT_k(T)$ (cf. [10], Chap. III, Prop. 3.6), we know that there is a suitable constant $\tilde{C}_{k,H}$ independent of $h$ such that

$$ \|\nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}]\|_{0,T} \le \bar{C}_{k,H} h^{k+1} \|\nabla \cdot \tilde{\mathbf{m}}\|_{k+1,T} \quad \forall T \in \mathcal{T}_h \setminus \mathcal{S}_{1,h}. \quad (161) $$

On the other hand, $\forall T \in \mathcal{S}_{1,h}$ we first write

$$ \left\{ \begin{array}{l} \|\nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}]\|_{0,T}^2 = (\nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}], \nabla \cdot \tilde{\mathbf{m}} - d_h(\tilde{\mathbf{m}}))_T \\ + (\nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}], d_h(\tilde{\mathbf{m}}) - \tilde{\Pi}_h \tilde{\mathbf{m}})_T. \end{array} \right. \quad (162) $$

Application of Young's inequality to (162), followed by the use of the second equation of (158), easily yields

$$ \left\{ \begin{array}{l} \|\nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}]\|_{0,T}^2 \le \tilde{C}_{k,d}^2 h^{2(k+1)} |\nabla \cdot \tilde{\mathbf{m}}|_{k+1,T^\Gamma}^2 + 2 \int_T \nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}] v \\ \text{with } v := d_h(\tilde{\mathbf{m}}) - \nabla \cdot \tilde{\Pi}_h \tilde{\mathbf{m}}. \end{array} \right. \quad (163) $$

Now, from the construction of $\tilde{\Pi}_h \tilde{\mathbf{m}}$, we readily establish that $\forall T \in \mathcal{S}_h$ and $\forall v \in P_k(T)$ it holds

$$ \int_T \nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}] v = \int_{e_T} [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}] \cdot \mathbf{n}_T v \, ds \quad \forall v \in P_k(T). \quad (164) $$

Let us set $\sigma(T) = 1$ if $T \in \mathcal{S}'_{1,h}$ and $\sigma(T) = -1$ if $T \in \mathcal{S}_{1,h} \setminus \mathcal{S}'_{1,h}$. Using the Divergence Theorem and taking into account that $\tilde{\mathbf{m}} \cdot \mathbf{n} = 0$ on $\Gamma_T$, after straightforward calculations starting from (164), we obtain $\forall v \in P_k(T^\Gamma)$,

$$ \left\{ \begin{array}{l} \int_T \nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}] v = I_{1,T}(\tilde{\mathbf{m}}, v) + I_{2,T}(\tilde{\mathbf{m}}, v) + I_{3,T}(\tilde{\mathbf{m}}, v) \\ \text{where} \\ I_{1,T}(\tilde{\mathbf{m}}, v) := -\sigma(T) \int_{\Gamma_T} [\tilde{\Pi}_h \tilde{\mathbf{m}}] \cdot \mathbf{n} \, v \, ds \\ I_{2,T}(\tilde{\mathbf{m}}, v) := \sigma(T) \int_{\Delta_T} [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}] \cdot \nabla v \\ I_{3,T}(\tilde{\mathbf{m}}, v) := \sigma(T) \int_{\Delta_T} \nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}] v. \end{array} \right. \quad (165) $$

Let us estimate $|I_{i,T}(\tilde{\mathbf{m}}, v)|$, $i = 1, 2, 3$.

First of all, recalling Lemma 4.3 and using the first equation of (158) with $j = k + 1$, we may write for a mesh-independent constant $C_{I1}$

$$ |I_{1,T}(\tilde{\mathbf{m}}, v)| \le C_{I1} h^{k+1} \|\tilde{\mathbf{m}}\|_{k+1+\varepsilon(1-\min[k,1]), T^\Gamma} \|v\|_{0,T}. \quad (166) $$

Next we turn our attention to $I_{2,T}(\tilde{\mathbf{m}}, v)$. Since this quantity vanishes for $k = 0$ $\forall v \in P_k(T)$, we handle it only for $k \ge 1$ as follows:

First we denote by $\tilde{\Pi}_T$ the counterpart in $\hat{T}^\Gamma$ of the interpolation operator $\tilde{\Pi}_h$ restricted to $T^\Gamma$ (cf. [10, 13]). In doing so, recalling the Piola transformation $\mathcal{P}_{T^\Gamma}$ from $[L^2(T^\Gamma)]^2$ onto $[L^2(\hat{T}^\Gamma)]^2$, we set $\tilde{\mathbf{m}} := \mathcal{P}_{T^\Gamma} \tilde{\mathbf{m}}$. Then denoting by $\rho_T$ the radius of the circle inscribed in $T$, referring to [13, 14], we observe that

$$ \|\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}\|_{0,\infty,\Delta_T} \le \left[ (2 - \sqrt{2}) \rho_T \right]^{-1} \|\tilde{\mathbf{m}} - \tilde{\Pi}_T \tilde{\mathbf{m}}\|_{0,\infty,\hat{T}^\Gamma}. \quad (167) $$
Now let $\hat{T}$ be the smallest isosceles right triangle that contains all $\hat{T}^\Gamma$ corresponding to $T \in S_{1,h}$ for all $T_h \in \mathcal{F}$ having the same right-angle vertex as all of them. It is not difficult to figure out that the diameter $\hat{h}$ of $\hat{T}$ is bounded above by $1 + \hat{c}h$, where $\hat{c}$ is a constant independent of $h$, and $\text{area}(\hat{T})$ is bounded below by $1/2$. Next, applying the inverse of $\mathcal{P}_{T^\Gamma}$ to $[L^2(\hat{T})]^2$ as well, we consider the extension $\hat{\mathbf{m}}$ of $\hat{\mathbf{m}}$ to $\hat{T} \setminus \hat{T}^\Gamma$ belonging to $[H^{k+1}(\hat{T})]^2$ given by $\hat{\mathbf{m}} = \mathcal{P}_{T^\Gamma} \tilde{\mathbf{m}}$.

Let us extend $\hat{\Pi}_T \hat{\mathbf{m}}$ to the whole $\hat{T}$. This extension in $\hat{T} \setminus \hat{T}^\Gamma$ equals $\hat{\Pi}_T \hat{\mathbf{m}}$ by construction. Indeed, the expression of $\hat{\Pi}_T \hat{\mathbf{m}}$ depends only on DOFs attached to $T^\Gamma$. Thus it is the same polynomial field as $\hat{\Pi}_T \hat{\mathbf{m}}$ in $\hat{T}^\Gamma$. It trivially follows that for a constant $C^\Gamma$ independent of $T^\Gamma$ it holds

$$ \| \hat{\Pi}_T \hat{\mathbf{m}} \|_{0,\infty, \hat{T}^\Gamma} \le C^\Gamma \| \hat{\mathbf{m}} \|_{0,\infty, \hat{T}^\Gamma} \le C^\Gamma \| \hat{\mathbf{m}} \|_{0,\infty, \hat{T}}. \quad (168) $$

Incidentally we observe that, by the Rellich-Kondratchev Theorem $H^{k+1}(\hat{T})$ is continuously embedded in $L^\infty(\hat{T})$ for $k \ge 1$. Moreover, since $\hat{\mathbf{m}} \equiv \hat{\Pi}_T \hat{\mathbf{m}}$ whenever $\hat{\mathbf{m}} \in P_k(\hat{T})$, by the Bramble-Hilbert lemma we can ensure the existence of a constant $\tilde{C}$ independent of $h$, such that

$$ \| \hat{\mathbf{m}} - \hat{\Pi}_T \hat{\mathbf{m}} \|_{0,\infty, \hat{T}^\Gamma} \le \| \hat{\mathbf{m}} - \hat{\Pi}_T \hat{\mathbf{m}} \|_{0,\infty, \hat{T}} \le \tilde{C} \| \hat{\mathbf{m}} \|_{k+1, \hat{T}}. \quad (169) $$

Now proceeding to the estimate of $I_{2,T}(\tilde{\mathbf{m}}, v)$, we first have

$$ |I_{2,T}(\tilde{\mathbf{m}}, v)| \le C_0 h_T^3 \|\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}\|_{0,\infty, \Delta_T} \|\nabla v\|_{0,\infty, T^\Gamma}. \quad (170) $$

Then, recalling (11), we obtain

$$ |I_{2,T}(\tilde{\mathbf{m}}, v)| \le C_0^+ C_0 h_T^2 \|\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}\|_{0,\infty, \Delta_T} \|\nabla v\|_{0,T}. \quad (171) $$

Further resorting to a standard inverse inequality, for a suitable constant $\tilde{C}_0$ independent of $T$ it holds

$$ |I_{2,T}(\tilde{\mathbf{m}}, v)| \le \bar{C}_0 h_T \|\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}\|_{0,\infty, \Delta_T} \|v\|_{0,T}. \quad (172) $$

Combining (167) and (169) and plugging the resulting inequality into (172) we come up with

$$ |I_{2,T}(\tilde{\mathbf{m}}, v)| \le \bar{C}_0 \tilde{C} h_T \left[ (2 - \sqrt{2}) \rho_T \right]^{-1} \|\hat{\mathbf{m}}\|_{k+1, \hat{T}} \|v\|_{0,T}. \quad (173) $$

Finally, let $\tilde{T} \subset \tilde{\Omega}$ be the reciprocal image of $\hat{T}$ under the same affine transformation from $T$ onto the standard master triangle $\hat{T}$. Using standard results for affine Piola transformations such as $\mathcal{P}_{T^\Gamma}$ (see e.g., [13, 14]), for a suitable mesh-independent constant $C_{I2}$, (173) easily leads to

$$ |I_{2,T}(\tilde{\mathbf{m}}, v)| \le C_{I2} h_T^{k+1} |\tilde{\mathbf{m}}|_{k+1, \tilde{T}} \|v\|_{0,T}. \quad (174) $$

Notice that several pairs of triangles $\tilde{T}$ overlap, but this does not affect the final result given below. As for $I_{3,T}(\tilde{\mathbf{m}}, v)$ we trivially have

$$ |I_{3,T}(\tilde{\mathbf{m}}, v)| \le [C_0 h_T^3]^{1/2} \|\nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}]\|_{0,T^\Gamma} \|v\|_{0,\infty, T^\Gamma}. \quad (175) $$

Recalling (11), for a suitable constant $\tilde{C}_{I3}$ the following estimate easily derives from (175):

$$ |I_{3,T}(\tilde{\mathbf{m}}, v)| \le \tilde{C}_{I3} h_T^{1/2} \|\nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}]\|_{0,T^\Gamma} \|v\|_{0,T}. \quad (176) $$
Now applying the triangle inequality and the second estimate of (158) to (176) we easily come up with a constant $C_{I3}$ such that

$$|I_{3,T}(\tilde{\mathbf{m}}, v)| \le C_{I3} h_T^{1/2} \left( h_T^{k+1} |\nabla \cdot \tilde{\mathbf{m}}|_{k+1,T^\Gamma} \|v\|_{0,T} + \|v\|_{0,T}^2 \right). \quad (177)$$

Plugging (166), (174), (177) into (165) and setting $\tilde{C}_I := \max[C_{I1}, C_{I2}, C_{I3}]$, we obtain $\forall T \in S_{1,h}$

$$\begin{aligned} \left\| \nabla \cdot (\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}) \right\|_{0,T}^2 &\le \tilde{C}_I h_T^{k+1} \left\{ \|\tilde{\mathbf{m}}\|_{k+1+\varepsilon(1-\min[k,1]),T^\Gamma} + |\tilde{\mathbf{m}}|_{k+1,\tilde{T}} + h_T^{1/2} |\nabla \cdot \tilde{\mathbf{m}}|_{k+1,T^\Gamma} \right\} \\ \|v\|_{0,T} + C_{I3} h_T^{1/2} \|v\|_{0,T}^2 &\le C_{I3} h_T^{1/2} \|v\|_{0,T}^2. \end{aligned} \quad (178)$$

Then, since $h_T < 1$, after straightforward manipulations we further establish the existence of a constant $\tilde{C}_I$ independent of $T$ such that

$$\begin{aligned} & (1 - 4C_{I3}h_T^{1/2}) \left\| \nabla \cdot (\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}) \right\|_{0,T}^2 \\ & \le \tilde{C}_I^2 h_T^{2k+2} \left\{ \|\tilde{\mathbf{m}}\|_{k+1+\varepsilon(1-\min[k,1]),T^\Gamma}^2 + |\tilde{\mathbf{m}}|_{k+1,\tilde{T}}^2 + |\nabla \cdot \tilde{\mathbf{m}}|_{k+1,T^\Gamma}^2 \right\}. \end{aligned} \quad (179)$$

Assuming that $h_T \le 1/(8C_{I3})^2$, we combine (179) and (161) to conclude that the interpolation error of $\nabla \cdot \tilde{\mathbf{m}}$ can be estimated by

$$\left\| \nabla \cdot [\tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}}] \right\|_{0,h} \le C_{k,div} h^{k+1} \left\{ \|\tilde{\mathbf{m}}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}} + \|\tilde{\mathbf{m}}\|_{k+1,\tilde{\Omega}} + |\nabla \cdot \tilde{\mathbf{m}}|_{k+1,\tilde{\Omega}} \right\}, \quad (180)$$

for a suitable constant $C_{k,div}$ independent of $h$.

Owing to (180) and (160), the result immediately follows with $C_H = [C_{k,div}^2 + C_{k,\tilde{\Omega}}^2]^{1/2}$. $\square$

# 7. ERROR ESTIMATES

In this section we endeavor to obtain error estimates for problem (7) based on the error bound (8). For the sake of clarity we distinguish the cases of convex and non convex domains.

## 7.1. The convex case

Here (8) can be applied with $\gamma_h = \gamma$ and $F_h(v) = (f, v)_h$. We consider separately the two possibilities, namely,

(I) $\Omega$ is convex and $length(\Gamma_0) > 0$

In this case we may take $(\tilde{\mathbf{p}}; \tilde{u}) = (\mathbf{p}; u)$.

The inf term in (8) can be bounded above by $|||(\mathbf{p}; u) - (\tilde{\Pi}_h \mathbf{p}; \tilde{\pi}_h u)|||_h$, where $\tilde{\pi}_h u$ is the standard Lagrange interpolate of $u$ for $k \ge 1$ and the mean value of $u$ in each triangle of the mesh for $k = 0$. We know that, as long as $u \in H^{k+1}(\Omega)$, there is a mesh-independent constant $C_\pi$ such that, for $k = 0, 1, \dots$ (cf. [13])

$$\|u - \tilde{\pi}_h u\|_{0,h} \le C_\pi h^{k+1} |u|_{k+1}. \quad (181)$$

On the other hand, the variational residual in the sup term of (8) is given by

$$\mathcal{R}_h(u, \mathbf{q}) := (\mathbf{p}, \mathbf{q})_h - (u, \nabla \cdot \mathbf{q})_h = -(u, \mathbf{q} \cdot \mathbf{n}_h)_{0,\Gamma_{0,h}} \quad \forall \mathbf{q} \in \mathbf{Q}_h. \quad (182)$$

In principle, $\mathcal{R}_h(u, \mathbf{q})$ does not vanish, since $u \equiv 0$ on $\Gamma_0$ but not necessarily on $\Gamma_{0,h}$. Nevertheless, we can give a rough estimate of this residual, in order to determine our method's real order in terms of $h$ if $length(\Gamma_0) > 0$.

**Proposition 7.1.** Let $length(\Gamma_0) > 0$ and $\mathcal{R}_h(u, \mathbf{q})$ be given by (182). If the first (resp. second) component $\mathbf{p}$ (resp. $u$) of the solution of (2) and (3) belongs to $[H^1(\Omega)]^2$ (resp. $H^2(\Omega)$) there exists a suitable mesh-independent constant $C_{\Gamma_0}$ such that for any $k \ge 0$ it holds

$$|\mathcal{R}_h(u, \mathbf{q})| \le C_{\Gamma_0} h^{3/2} \|u\|_2 [|\mathbf{q}|]_h \quad \forall \mathbf{q} \in \mathbf{Q}_h. \quad (183)$$
*Proof.* First of all we have

$$ (u, \mathbf{q} \cdot \mathbf{n}_h)_{e_T} \leq \|u\|_{0, e_T} \|\mathbf{q} \cdot \mathbf{n}_h\|_{0, e_T} \quad \forall T \in \mathcal{S}_{0, h}. \qquad (184) $$

By a straightforward manipulation already performed several times in this work, we come up with a constant $C_{0I}$ independent of $T$ such that

$$ \|\mathbf{q} \cdot \mathbf{n}_h\|_{0, e_T} \le l_T^{1/2} \|\mathbf{q}\|_{0, \infty, T} \le C_{0I} h_T^{-1/2} \|\mathbf{q}\|_{0, T}. \qquad (185) $$

Moreover, since $u$ vanishes at both ends of $e_T$, we know that, for a certain constant $C_{e0}$, it holds

$$ \|u\|_{0, e_T} \le C_{e0} h_T^j \|\partial^j u / \partial x^j\|_{0, e_T} \quad \text{for } 1 \le j \le 2. \qquad (186) $$

Hereafter we take $j = 1$. Recalling the local coordinate system $(A_T; x, y)$ attached to $e_T$ and the function $\phi$ such that $M = (x, 0) \in e_T$ corresponds to $N = (x, \phi(x)) \in \Gamma_T$. We have

$$ (\partial u / \partial x)(x, 0) = (\partial u / \partial x)(x, \phi(x)) - \int_M^N \partial^2 u / \partial x \partial y \, dy. \qquad (187) $$

$\mathbf{t}$ being the tangent unit vector along $\Gamma$ we introduce the direct reference frame $\mathcal{R} := (N; \mathbf{t}, \mathbf{n})$. Then recalling the angle $\varphi$ between $\mathbf{n}_h$ and $\mathbf{n}$, let the unit vector along $e_T$ be $\mathbf{l}_T$. We have $\mathbf{l}_T = \cos \varphi \mathbf{t} - \sin \varphi \mathbf{n}$ and hence we have $\partial u / \partial x = \nabla u |_{\mathcal{R}} \cdot (\cos \varphi; -\sin \varphi)$. We know that $|\sin \varphi| \le C_T h_T$ owing to (9) and also that all the tangential derivatives of $u$ of any order vanish on $\Gamma_0$. This easily yields

$$ |(\partial u / \partial x)(x, \phi(x))| \le C_T h_T |(\partial u / \partial n)(x, \phi(x))| \quad \forall x \in [0, l_T]. \qquad (188) $$

Plugging (188) into (187) we obtain

$$ |(\partial u / \partial x)(x, 0)| \le C_T h_T |(\partial u / \partial n)(x, \phi(x))| + C_0^{1/2} h_T \left[ \int_M^N |\partial^2 u / \partial x \partial y|^2 \, dy \right]^{1/2} \quad \forall x \in [0, l_T]. \qquad (189) $$

It is easy to see that $\cos \varphi \, ds = dx$ where $s$ is the curvilinear abscissa along $\Gamma_T$, that is, $ds$ is just a small perturbation of $dx$. Thus plugging in turn (189) into (186), it trivially holds for a certain constant $C_{e1}$ independent of $T$:

$$ \|u\|_{0, e_T}^2 \le C_{e1}^2 h_T^4 \left[ \|\partial u / \partial n\|_{0, \Gamma_T}^2 + |u|_{2, \Delta_T}^2 \right]. \qquad (190) $$

Finally recalling (184), (185) and sweeping the triangles $T \in \mathcal{S}_{0, h}$, straightforward manipulations yield

$$ (u, \mathbf{q} \cdot \mathbf{n}_h)_{0, \Gamma_{0,h}} \le C_{0I} C_{e1} h^{3/2} \{ \|\partial u / \partial n\|_{0, \Gamma_0} + |u|_2 \} \|\mathbf{q}\|_h. \qquad (191) $$

Finally, thanks to the Trace Theorem for $H^2(\Omega)$ (cf. [1]), (183) derives from (191). $\square$

As a consequence of Proposition 7.1 we trivially have

**Theorem 7.2.** Assume that length($\Gamma_0$) > 0, $\Omega$ is smooth enough, $h$ is sufficiently small and $k \ge 0$. If for a small $\varepsilon > 0$ $\tilde{\mathbf{p}} \in [H^{k+1+\varepsilon(1-\min[1,k])}(\tilde{\Omega})]^2$, $\nabla \cdot \mathbf{p} \in H^{k+1}(\Omega)$ and $u \in H^{\max[2,k+1]}(\Omega)$, the following error estimate holds with a mesh-independent constant $C_{01}$, for the approximation of the solution of (2) and (3) by the solution of (7).

$$ |||(\mathbf{p}; u) - (\mathbf{p}_h; u_h)|||_h \le C_{01} \left\{ h^{3/2} \|u\|_2 + h^{k+1} \left[ \|\tilde{\mathbf{p}}\|_{k+1+\varepsilon(1-\min[1,k]), \tilde{\Omega}} + |\nabla \cdot \mathbf{p}|_{k+1} + |u|_{k+1} \right] \right\}. \qquad (192) $$

**Remark 7.3.** According to the estimate (192), in case length($\Gamma_0$) > 0, the order of our method is optimal (up to a very small number $\varepsilon$) only for $k = 0$. But there should be no grievance about this fact, since in general the solution $u$ of the mixed Poisson problem does not even belong to $H^2(\Omega)$, but rather to $H^{3/2+\varepsilon}(\Omega)$ [22]. Therefore $\mathbf{p} \in [H^{1/2+\varepsilon}(\Omega)]^2$ in general, and hence the best order we can hope for here is a little more than 1/2. Nevertheless there is a rather simple way out, if the exact solution happens to be more regular, as we will see in Section 7.3. Incidentally, we will also see later on that, by requiring more regularity from $u$, our method’s optimal (second) order can be recovered for $k = 1$, in the case of mixed boundary conditions.
(II) $\Omega$ is convex and length($\Gamma_0$) = 0

Here it is advisable to take $(\tilde{\mathbf{p}}; \tilde{u}) = (\mathbf{p}; \bar{u})$ where $\bar{u}$ differs from $u$ by a suitable additive constant. Indeed, since $u \in L_0^2(\Omega)$ in this case, a natural way to define the interpolate of $u$ in $V_h$ is as follows: We start with the operator $\tilde{\pi}_h$ from $V^k(\Omega) := H^{k+1}(\Omega) \cap V$ into $L_0^2(\Omega_h)$ already defined in I). Noticing that $\tilde{\pi}_h u$ does not necessarily belong to $L_0^2(\Omega_h)$, $\tilde{\pi}_h u$ may not belong to $V_h$. Thus we set $\tilde{\pi}_h u := \tilde{\pi}_h u + C_h \in V_h$ with $C_h := - \int_{\Omega_h} \tilde{\pi}_h u/\text{area}(\Omega_h)$. While on the one hand the estimate of the interpolation error in $V_h$ of $u$ is not optimal, we note that the estimate $\|u - \tilde{\pi}_h u\|_{0,h} \le C_\pi h^{k+1} |u|_{k+1}$ holds if $u \in V^k(\Omega)$. Now we define $\bar{u} = u + C_h$, which leads to the desired optimal estimate of the best approximation in $\mathbf{P}_h \times V_h$ of the pair $(\mathbf{p}; \bar{u})$ instead of $(\mathbf{p}; u)$, as seen below.

It is noticeable that, since $\int_{\Omega} \nabla \cdot \mathbf{q} = 0 \quad \forall \mathbf{q} \in \mathbf{Q}$, the replacement of $u$ by $\bar{u}$ does not change the essence of problem (2) and (3); it is just a way like any other to fix the additive constant up to which $u$ is defined.

Summarizing, the error bound (8) applies to the pair $(\tilde{\mathbf{p}}; \tilde{u}) = (\mathbf{p}; \bar{u})$. It turns out that

$$ \inf_{w \in V_h} \|\bar{u} - w\|_{0,h} \le \|\bar{u} - \tilde{\pi}_h u\|_{0,h} = \|u - \tilde{\pi}_h u\|_{0,h} \le C_\pi h^{k+1} |u|_{k+1}. \quad (193) $$

Thus recalling (159), provided $\tilde{\mathbf{p}} \in [H^{k+1+\varepsilon(1-\min[k,1])}(\tilde{\Omega})]^2$ and $\nabla \cdot \mathbf{p} \in H^{k+1}(\Omega)$, setting $C_I := C_H + C_\pi$, we come up with,

$$ \inf_{(\mathbf{r};w) \in \mathbf{P}_h \times V_h} |||(\mathbf{p}; \bar{u}) - (\mathbf{r}; w)|||_h \le C_I h^{k+1} \left\{ \|\tilde{\mathbf{p}}\|_{k+1+\varepsilon(1-\min[k,1]), \tilde{\Omega}} + |\nabla \cdot \mathbf{p}|_{k+1} + |u|_{k+1} \right\}. \quad (194) $$

As for the $\sup$ term in (8), we can assert that it vanishes identically if $\text{length}(\Gamma_0) = 0$. Indeed, $\nabla \cdot \mathbf{p} = -f$ in $\Omega_h$ and on the other hand, as long as $u \in H^{k+1}(\Omega)$,

$$ (\mathbf{p}, \mathbf{q})_h - (\bar{u}, \nabla \cdot \mathbf{q})_h = (\mathbf{p} - \nabla \bar{u}, \mathbf{q})_h = 0 \quad \forall \mathbf{q} \in \mathbf{Q}_h, \quad (195) $$

since $\mathbf{q} \cdot \mathbf{n}_h$ vanishes everywhere in $\Gamma_h$ in this case.

This leads to the following result.

**Theorem 7.4.** Assume that $\text{length}(\Gamma_0) = 0$, $\Omega$ is smooth enough and $h$ is sufficiently small. If $\tilde{\mathbf{p}} \in [H^{k+1+\varepsilon(1-\min[k,1])}(\tilde{\Omega})]^2$, $\nabla \cdot \mathbf{p} \in H^{k+1}(\Omega)$ and $u \in H^{k+1}(\Omega)$, the following error estimate holds for the approximation of the solution of (2) and (3) by the solution of (7).

$$ \left\{ \begin{array}{l} |||(\mathbf{p}; \bar{u}) - (\mathbf{p}_h; u_h)|||_h \le \frac{C_I}{\gamma} h^{k+1} \left\{ \|\tilde{\mathbf{p}}\|_{k+1+\varepsilon(1-\min[k,1]), \tilde{\Omega}} + |\nabla \cdot \mathbf{p}|_{k+1} + |u|_{k+1} \right\} \\ \text{where } u - \bar{u} = \text{constant, for } k = 0, 1, 2, \dots \end{array} \right. \quad (196) $$

## 7.2. The non convex case

Here we can apply (8) with $\gamma_h = \gamma'$, (cf. (155)) but $F_h(v)$ cannot be $(f, v)_h$. Besides this one, new issues must be overcome in the non convex case as far as error estimates are concerned.

First of all, whenever $\Omega$ is not convex, the error cannot be measured in norms associated with the set $\Omega_h$, since the exact solution is not defined in $\Omega_h \setminus \tilde{\Omega} \neq \emptyset$. Akin to [27] the error estimate is now expressed for $|||(\mathbf{p}; u) - (\mathbf{p}_h; u_h)|||'_h$ which is not greater than $|||(\tilde{\mathbf{p}}; \tilde{u}) - (\mathbf{p}_h; u_h)|||_h$ for a suitable pair $(\tilde{\mathbf{p}}; \tilde{u})$ defined hereafter. This is in turn bounded above as in (8).

As for functional $F_h$, it must be defined in such a way that no unknown values of an extension $\tilde{f}$ of $f$ to $\Omega_h \setminus \Omega$ come into play. Notice that the replacement of $(\tilde{f}, v)_h$ by $F_h(v)$ on the right hand side of (7) is an old trick considered in [15] to handle curved non convex domains.

Incidentally we refer to Section 5.3 for some definitions and concepts evoked hereafter concerning extensions to $\tilde{\Omega}$, i.e., a smooth open set sufficiently large to strictly contain not only $\Omega$, but also the sets $\tilde{\Omega}_h$ for all $\mathcal{T}_h \in \mathcal{F}'$.

Actually the pair $(\tilde{\mathbf{p}}; \tilde{u})$ extending $(\mathbf{p}; u)$ to $\tilde{\Omega}$ must have the required regularity for method's order to remain the same as in the convex case. Akin to [27], the construction of the extension $\tilde{u}$ of $u$ which incidentally induces regular extensions $\tilde{\mathbf{p}} := \nabla \tilde{u}$ and $\tilde{f} := -\Delta \tilde{u}$ of $\mathbf{p}$ and $f$, is the one advocated in [33]. Of course, the
explicit knowledge of neither $(\tilde{\mathbf{p}}; \tilde{u})$ nor $\tilde{f}$ is necessary. The point here is that all these quantities do exist, and can thus be employed in error estimates for the non convex case, which we next complete.

Notice that the extension technique allows for establishing an estimate for the interpolation error in the non convex case analogous to (159). More precisely if $\tilde{\mathbf{m}} \in \mathcal{M}^k(\tilde{\Omega})$ we have

$$ \left[ \left| \tilde{\mathbf{m}} - \tilde{\Pi}_h \tilde{\mathbf{m}} \right| \right]_h \le C_H h^{k+1} \left[ \| \tilde{\mathbf{m}} \|_{k+1+\varepsilon(1-\min[k,1]), \tilde{\Omega}}^2 + |\nabla \cdot \tilde{\mathbf{m}}|_{k+1, \tilde{\Omega}}^2 \right]^{1/2} \quad (197) $$

Furthermore, we may exploit the extension $\tilde{f}$ in the framework of the approximate right hand side $F_h(v)$ of problem (7). Among many other possibilities, for the sake of simplicity we consider here a bilinear form $F_h$ defined as follows. First of all, assuming that $\tilde{f} \in H^{k+1+\varepsilon(1-\min[k,1])}(\tilde{\Omega}_h)$, let $f_h$ be the Lagrange $P_k(T)$- interpolate of $\tilde{f}$ at the set of $(k+2)(k+1)/2$ points in $T$, $\forall T \in \mathcal{T}_h$, which form the standard set of lattice points of $T_\chi$, where $T_\chi$ is the reduction of $T$ by homothety with ratio $\chi < 1$ centered at the centroid of $T$. $\chi$ is chosen in such a way that $T_\chi \subset \tilde{\Omega} \ \forall T \in \mathcal{S}_h$ for all $\mathcal{T}_h$ in $\mathcal{F}'$.

Then we set

$$ F_h(v) := (f_h, v)_h \quad \forall v \in V_h. \quad (198) $$

Since by standard results [13] there is a constant $C_{\chi,k}$ independent of $h$ such that

$$ \| \tilde{f} - f_h \|_{0,h} \le C_{\chi,k} h^{k+1} \| \tilde{f} \|_{k+1+\varepsilon(1-\min[k,1]), \Omega_h}, $$

we obtain from the Cauchy-Schwarz inequality, together with (198)

$$ \left| (f, v)_h - F_h(v) \right| \le C_{\chi,k} h^{k+1} \left\| \tilde{f} \right\|_{k+1+\varepsilon(1-\min[k,1]), \Omega_h} \|v\|_{0,h} \quad \forall v \in V_h. \quad (199) $$

Henceforth we definitively take $\tilde{f} = -\Delta \tilde{u}$ in $\tilde{\Omega} \setminus \Omega$ as the extension of $f$ with the required smoothness.

In view of all the above considerations, we are now ready to establish error estimates for problem (7) in the case of non convex domains. Akin to the previous subsection we study separately the two possible cases, namely,

### (III) $\tilde{\Omega}$ is not convex and $\text{length}(\Gamma_0) > 0$

We next focus on the estimation of the contribution to the approximation error of the variational residual incorporating the functional $F_h(v)$, namely, $|((\tilde{\mathbf{p}}; \tilde{u}), (\mathbf{q}; v)) - F_h(v)|$. With this aim we first write

$$ |c_h((\tilde{\mathbf{p}}; \tilde{u}), (\mathbf{q}; v)) - F_h(v)| \le |c_h((\tilde{\mathbf{p}}; \tilde{u}), (\mathbf{q}; v)) - (\tilde{f}, v)_h| + |(\tilde{f}, v)_h - F_h(v)|. \quad (200) $$

Now we observe that, either by construction, or due to (1), it holds $-\nabla \cdot \tilde{\mathbf{p}} = \tilde{f}$ and $\tilde{\mathbf{p}} = \nabla \tilde{u}$ almost everywhere in $\tilde{\Omega}$. Thus assuming that $\tilde{u} \in H^2(\tilde{\Omega})$, the first term on the right hand side of (200) is $\mathcal{R}_h(\tilde{u}, \mathbf{q}) := (\tilde{u}, \mathbf{q} \cdot \mathbf{n}_h)_{0,\Gamma_{0,h}}$. It is easy to see that $\mathcal{R}_h(\tilde{u}, \mathbf{q})$ can be estimated in the same manner as in (183), that is,

$$ |\mathcal{R}_h(\tilde{u}, \mathbf{q})| \le C_{\Gamma_0} h^{3/2} \|\tilde{u}\|_{2,\tilde{\Omega}_h} |[\mathbf{q}]|_h \quad \forall \mathbf{q} \in \mathbf{Q}_h. \quad (201) $$

Combining (199) and (201) and setting $\tilde{C}_r := \max[C_{\Gamma_0}, C_{\chi,k}]$, we immediately come up with

$$ \sup_{(\mathbf{q}; v) \in \mathbf{Q}_h \times V_h / |||(\mathbf{q}; v)|||_h = 1} |c_h((\tilde{\mathbf{p}}; \tilde{u}), (\mathbf{q}; v)) - F_h(v)| \le \tilde{C}_r \left[ h^{3/2} \|\tilde{u}\|_{2,\tilde{\Omega}} + h^{k+1} \|\tilde{f}\|_{k+1+\varepsilon(1-\min[k,1]), \tilde{\Omega}} \right]. \quad (202) $$

Now taking $\tilde{\mathbf{m}} = \tilde{\mathbf{p}}$ in (197) and letting $\tilde{\pi}_h \tilde{u}$ be the $V_h$-interpolate of $\tilde{u}$ for $k > 0$ and the mean value of $\tilde{u}$ in every triangle of the mesh for $k = 0$, we can assert that there exists a constant $\tilde{C}_I$ independent of $h$ such that the following estimate holds:

$$ \inf_{(\mathbf{r}; w) \in \mathbf{P}_h \times V_h} |||(\tilde{\mathbf{p}}; \tilde{u}) - (\mathbf{r}; w)|||_h \le |||(\tilde{\mathbf{p}}; \tilde{u}) - (\tilde{\Pi}_h \tilde{\mathbf{p}}; \tilde{\pi}_h \tilde{u})|||_h $$
$$ \leq \tilde{C}_I h^{k+1} \left\{ \|\tilde{\mathbf{p}}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}} + |\nabla \cdot \tilde{\mathbf{p}}|_{k+1,\tilde{\Omega}} + |\tilde{u}|_{k+1,\tilde{\Omega}} \right\}. \quad (203) $$

Finally, plugging (202) and (203) into (8), we obtain an error estimate for the non convex case with mixed boundary conditions. More precisely, we have

**Theorem 7.5.** Assume that $\Omega$ is not convex and length($\Gamma_0$) > 0. If $\mathbf{p} \in [H^{k+1+\varepsilon(1-\min[k,1])}(\Omega)]^2$, $u \in H^{\max[2,k+1]}(\Omega)$ and $\nabla \cdot \mathbf{p} \in H^{k+1}(\Omega)$, then

$$ \begin{aligned} |||(\mathbf{p}; u) - (\mathbf{p}_h; u_h)|||_h' &\leq \max[\tilde{C}_I, \tilde{C}_r] \left\{ h^{3/2} \|\tilde{u}\|_{2,\tilde{\Omega}} \right. \\ & \qquad \left. + h^{k+1} \left[ \|\tilde{f}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}} + \|\tilde{\mathbf{p}}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}} + |\nabla \cdot \tilde{\mathbf{p}}|_{k+1,\tilde{\Omega}} + |\tilde{u}|_{k+1,\tilde{\Omega}} \right] \right\}. \end{aligned} \quad (204) $$

From (204) we infer that, in case mixed conditions are prescribed on the boundary of a non convex domain, our method in original version enjoys an optimal order of convergence only for $k = 0$, even if at the price of the stronger regularity assumption $u \in H^2(\Omega)$.

We next explain why an optimal order can also be expected of our method for $k = 1$, in case $u \in H^3(\Omega)$. With this aim we first note that

$$ |\mathcal{R}_h(\tilde{u}, \mathbf{q})| \leq \|\tilde{u}\|_{1/2, \Gamma_{0,h}} \|\mathbf{q} \cdot \mathbf{n}_h\|_{-1/2, \Gamma_{0,h}}, \quad (205) $$

where $\| \cdot \|_{-1/2, D}$ represents the norm of the dual space of $H^{1/2}(D)$, for any $D \subset \tilde{\Omega}$. According to [21] there is a constant $C_{\Gamma_{0,h}}$ depending on $\Omega_h$ such that

$$ \|\mathbf{q} \cdot \mathbf{n}_h\|_{-1/2, \Gamma_{0,h}} \leq C_{\Gamma_{0,h}} [\|\mathbf{q}\|_h]. \quad (206) $$

On the other hand we have

$$ \|\tilde{u}\|_{1/2, \Gamma_{0,h}} = \left[ \sum_{T \in S_{0,h}} \|\tilde{u}\|_{1/2, e_T}^2 \right]^{1/2} . \quad (207) $$

It turns out that $\tilde{u}$ vanishes at the end points of $e_T$. Therefore, by standard approximation results (cf. [32]) there is a mesh-independent constant $C_u$ such that

$$ \|\tilde{u}\|_{1/2, e_T} \leq C_u l_T^2 |\tilde{u}|_{5/2, e_T}. \quad (208) $$

Plugging (208) into (207) and combining the resulting expression with (206), after straightforward calculations we conclude from (205) that

$$ |\mathcal{R}_h(\tilde{u}, \mathbf{q})| \leq C_u C_{\Gamma_{0,h}} h^2 |\tilde{u}|_{5/2, \Gamma_{0,h}} [\|\mathbf{q}\|_h]. \quad (209) $$

Now, by the Trace Theorem, there exists another constant $C_{\Omega_h}$ depending only on $\Omega_h$ such that

$$ |\tilde{u}|_{5/2, \Gamma_{0,h}} \leq C_{\Omega_h} \|\tilde{u}\|_{3, \Omega_h}. \quad (210) $$

Plugging (210) into (209) we conclude that the residual $\mathcal{R}_h(\tilde{u}, \mathbf{q})$ is likely to be an $\mathcal{O}(h^2)$-term. More precisely, there is a constant $C_{\tilde{\Omega}_h}$ such that

$$ |\mathcal{R}_h(\tilde{u}, \mathbf{q})| \leq C_{\tilde{\Omega}_h} h^2 \|\tilde{u}\|_{3, \tilde{\Omega}} [\|\mathbf{q}\|_h]. \quad (211) $$

Due to the closeness of $\Gamma_0$ and $\Gamma_{0,h}$, it can be strongly conjectured that $C_{\tilde{\Omega}_h}$ is mesh-independent. However we refrain from elaborating any further on this issue for the sake of brevity. Actually, in Appendix II of [4] we show through a numerical example that, in the case of mixed boundary conditions too, our method is indeed of the second order for $k = 1$, as long as $u$ is sufficiently smooth.
(IV) $\Omega$ is not convex and $length(\Gamma_0) = 0$

First of all we note that here, while $\tilde{\mathbf{p}}$ in (8) is the extension of $\mathbf{p}$ to $\tilde{\Omega}$, akin to the convex case, $\tilde{u}$ is not the aforementioned extension of $u$ to the same set, but rather a function $\bar{u}$ which differs from it by a suitable additive constant. Indeed, recalling the operator $\tilde{\pi}_h$ defined in the previous subsection, first we introduce $\bar{u} = \tilde{u} + C_h$ in $\tilde{\Omega}$, where $C_h := -\int_{\Omega_h} \tilde{\pi}_h \tilde{u}/\text{area}(\Omega_h)$. Then we note that the estimate $\|\tilde{u} - \tilde{\pi}_h \tilde{u}\|_{0,h} \le C_\pi h^{k+1}|\tilde{u}|_{k+1,\tilde{\Omega}}$ holds if $u \in H^{k+1}(\Omega)$, but not necessarily for $u \in V^k(\Omega)$. Therefore, defining in $\Omega_h$ $\bar{\pi}_h \tilde{u} := \tilde{\pi}_h \tilde{u} + C_h$ we have $\bar{\pi}_h \tilde{u} \in V_h$ and thus

$$ \left\{ \begin{array}{ll} \|\bar{u} - \bar{\pi}_h \tilde{u}\|_{0,h} \le C_\pi h^{k+1} |\tilde{u}|_{k+1,\tilde{\Omega}} & \text{where} \\ \bar{u} = \tilde{u} + C_h \text{ and } \bar{\pi}_h \tilde{u} := \tilde{\pi}_h \tilde{u} + C_h & \\ \text{with } C_h := -\int_{\Omega_h} \tilde{\pi}_h \tilde{u}/\text{area}(\Omega_h). & \end{array} \right. \qquad (212) $$

Equation (212) can be used to estimate the inf-term in (8), that is,

$$ \inf_{w \in V_h} \|\bar{u} - w\|_{0,h} \le \|\bar{u} - \bar{\pi}_h \tilde{u}\|_{0,h} = \|\tilde{u} - \tilde{\pi}_h \tilde{u}\|_{0,h} \le C_\pi h^{k+1} |\tilde{u}|_{k+1, \tilde{\Omega}}. \qquad (213) $$

**Remark 7.6.** $\bar{u}$ can be replaced by the regular extension $\tilde{u}$ of $u$ in the sup-term of (8) since $(C_h, \nabla \cdot \mathbf{q})_h = 0$ $\forall \mathbf{q} \in \mathbf{Q}_h$ if $length(\Gamma_0) = 0$.

Recalling the treatment of the case $length(\Gamma_0) = 0$ for convex domains, an error estimate like (204) holds for non convex domains without the term in $h^{3/2}$, though for a pair $(\mathbf{p}; \bar{u})$ instead of $(\mathbf{p}; u)$, where $\bar{u}$ is a function that differs from $u$ by a suitable additive constant. Indeed now sole $|(\tilde{f}, v)_h - F_h(v)|$ appears in the sup-term of (8). This leads to

**Theorem 7.7.** Assume that $\Omega$ is not convex and $length(\Gamma_0) = 0$. If $\mathbf{p} \in [H^{k+1+\varepsilon(1-\min[k,1])}(\Omega)]^2$, $u \in H^{k+1}(\Omega)$ and $\nabla \cdot \mathbf{p} \in H^{k+1}(\Omega)$, then setting $\bar{u} = u + C_h$ we have

$$ \|\|(\mathbf{p}; \bar{u}) - (\mathbf{p}_h; u_h)\|\|'_h \le \max[\tilde{C}_I, C_{\chi,h}] h^{k+1} \\ \times \left[ \|\tilde{f}\|_{k+1+\varepsilon(1-\min[k,1]), \tilde{\Omega}} + \|\tilde{\mathbf{p}}\|_{k+1+\varepsilon(1-\min[k,1]), \tilde{\Omega}} + |\nabla \cdot \tilde{\mathbf{p}}|_{k+1, \tilde{\Omega}} + |\tilde{u}|_{k+1, \tilde{\Omega}} \right]. \quad (214) $$

Theorem 7.7 indicates that, in case pure Neumann conditions are prescribed on the boundary of a non convex domain, our method enjoys an optimal order of convergence for every $k \ge 0$, provided the exact solution of (2) and (3) is sufficiently smooth.

**Remark 7.8.** If a convergence check is being carried out for a test-problem having a known solution $u$ sufficiently smooth, we can certainly exhibit $\tilde{u}$ in $\tilde{\Omega}_h$ for every triangulation of interest. But in this case it is simpler to define $\bar{u}$ by $u - u(Q)$ where $Q$ is a well chosen point in $\tilde{\Omega}_h$ at which we prescribe $u_h(Q) = 0$ for all the triangulations of interest. In this manner the error $\|\bar{u} - u_h\|'_{0,h}$ – or yet, more conveniently, $\|[\tilde{u} - u(Q)] - u_h\|_{0,h}$ –, is the one to be computed for every $h$ under consideration, for it will decrease at the best rate.

### 7.3. Recovering optimal orders for mixed boundary conditions with arbitrary $k$

As pointed out in Remark 7.3 of Section 7.1, the error estimate (192) applying to the case of mixed boundary conditions and convex domains is somewhat optimistic even for $k = 0$. However in the non convex case this is no longer true if the intersection between $\bar{\Gamma}_0$ and $\bar{\Gamma}_1$ is empty. This is because in this case the exact solution of (1) can be arbitrarily smooth, as long as $\Gamma$ and $f$ are sufficiently regular. As we saw in Section 7.1, the problem stems from the term in the variational residual involving an integral along $\Gamma_{0,h}$, which is an $O(h^{3/2})$ whatever $k$. Nevertheless this issue can be settled for $k > 0$ if every (straight) triangle in $T \in S_{0,h}$ is replaced by a triangle $\check{T}$ having the same straight edges as $T$ except $e_T$, which is replaced by a curvilinear edge $\check{e}_T$ having the same
end-points as $e_T$. The idea here is not to use a parametric representation, since the master element $\tilde{T}$ associated with $\tilde{T}$ is also a triangle with a curvilinear edge. Actually $\tilde{T}$ results from the same affine mapping as the one from $T$ onto the standard master element $\hat{T}$. Therefore the integrations can be easily performed, either in $\tilde{T}$ itself or in $\hat{T}$. The point here is that the curvilinear edge $\check{e}_T$ is generated as the polynomial arc of degree $k+1$ which interpolates $\Gamma_0$ at $k+2$ points including the end points of $e_T$. In this manner, the residual term involving the integral of $u\mathbf{q} \cdot \mathbf{n}_h$ on $\Gamma_{0,h}$, with $\mathbf{q} \in \mathbf{Q}_h$, becomes

$$ \check{\mathcal{R}}_h(\tilde{u}, \mathbf{q}) := \sum_{T \in \mathcal{S}_{0,h}} (\tilde{u}, \mathbf{q} \cdot \check{\mathbf{n}}_T)_{\check{e}_T}, \qquad (215) $$

$\mathbf{q}$ being extended to a non empty $\tilde{T} \setminus T$, where $\check{\mathbf{n}}_T$ is the outer normal vector to $\check{e}_T$.

Likewise the convex case (for $k \ge 1$), we assume that $u \in H^3(\Omega)$, or yet that $\tilde{u} \in H^3(\tilde{\Omega})$.

Owing to the construction of $\check{e}_T$, $\check{\mathcal{R}}_h(\tilde{u}, \mathbf{q})$ can be estimated as an $O(h^{k+3/2})$ provided $\tilde{u} \in H^3(\tilde{\Omega})$. Indeed, akin to (184) and (185), we first establish the existence of a constant $\check{C}_{0I}$ independent of $T$ such that

$$ (\tilde{u}, \mathbf{q} \cdot \check{\mathbf{n}}_h)_{\check{e}_T} \le \check{C}_{0I} h_T^{-1/2} \|\tilde{u}\|_{0, \check{e}_T} \|\mathbf{q}\|_{0, T} \quad \forall T \in \mathcal{S}_{0,h}. \qquad (216) $$

Therefore, recalling (215) we have

$$ \check{\mathcal{R}}_h(\tilde{u}, \mathbf{q}) \le \check{C}_{0I} \left[ \sum_{T \in \mathcal{S}_{0,h}} h_T^{-1} \|\tilde{u}\|_{0, \check{e}_T}^2 \right]^{1/2} \|\mathbf{q}\|_h. \qquad (217) $$

Now let $P$ be the point of $\check{e}_T$ located on the perpendicular to $e_T$ passing through a given point $M \in e_T$. Referring to Figure 1 and recalling the local coordinate system $(A_T; x, y)$ attached to $e_T$, since $\tilde{u}(N) = 0$, whatever $M \in e_T$ we have

$$ \tilde{u}(P) = \int_N^P \frac{\partial \tilde{u}}{\partial y} dy \quad \forall P \in \check{e}_T. \qquad (218) $$

Let $\phi_h(x)$ be the function that uniquely represents $\check{e}_T$ in such a frame. The curvilinear abscissa $s_h$ along $\check{e}_T$ satisfies $ds_h = \sqrt{1 + [\phi'_h(x)]^2} dx$. On the other hand, by basic interpolation theory in one-dimensional space (see e.g., [7,24]), there exists a constant $C_L$ independent of $l_T$ such that $\|\phi'_h\|_{0,\infty,e_T} \le C_L \|\phi'\|_{0,\infty,e_T} \le C_\Gamma C_L h_T$. Hence from (218) we easily establish the existence of a constant $\check{C}_L$ independent of $T$ fulfilling

$$ \int_{\check{e}_T} \tilde{u}^2(P) ds_h = \int_{\check{e}_T} \left[ \int_N^P \frac{\partial \tilde{u}}{\partial y} dy \right]^2 ds_h \le \check{C}_L \int_0^{l_T} \left[ \int_N^P \frac{\partial \tilde{u}}{\partial y} dy \right]^2 dx. \qquad (219) $$

Taking an obvious upper bound of the right hand side of (219), we readily obtain $\forall T \in \mathcal{S}_{0,h}$

$$ h_T^{-1} \|\tilde{u}\|_{0, \check{e}_T}^2 \le \check{C}_L \max_{x \in [0, l_T]} |\phi_h(x) - \phi(x)|^2 \|\tilde{u}\|_{1, \infty, \check{\Omega}}^2. \qquad (220) $$

Resorting again to the interpolation theory in one-dimensional space, since a sufficient smoothness of $\Omega$ is assumed, we can assert that there are two constants $C_{L,k}$ and $\check{C}_k$, both independent of $e_T$, such that

$$ \max_{x \in [0, l_T]} |\phi_h(x) - \phi(x)| \le C_{L,k} l_T^{k+2} |\phi|_{k+2, \infty, e_T} \le \check{C}_k l_T^{k+2}. \qquad (221) $$

Plugging (221) into (220) and taking into account the continuous embedding of $H^3(\tilde{\Omega})$ into $W^{1,\infty}(\tilde{\Omega})$, summing up over $\mathcal{S}_{0,h}$ we readily come up with a mesh-independent constant $C_{\tilde{\Omega}}$ satisfying

$$ \sum_{T \in \mathcal{S}_{0,h}} h_T^{-1} \|\tilde{u}\|_{0, \check{e}_T}^2 \le \check{C}_L \check{C}_k^2 h^{2k+3} \|\tilde{u}\|_{1, \infty, \check{\Omega}}^2 \sum_{T \in \mathcal{S}_{0,h}} l_T \le C_{\tilde{\Omega}}^2 \text{length}(\Gamma_0) h^{2k+3} \|\tilde{u}\|_{3, \check{\Omega}}^2. \qquad (222) $$
Finally, combining (222) and (217), we infer the existence of a constant $\check{C}_{\Gamma_0}$ depending only on $\Gamma_0$ and $\check{\Omega}$ such that the residual term $\check{\mathcal{R}}_h(\tilde{u}, \mathbf{q})$ fulfills

$$ \left| \check{\mathcal{R}}_h(\tilde{u}, \mathbf{q}) \right| \le \check{C}_{\Gamma_0} h^{k+3/2} \|\tilde{u}\|_{3, \check{\Omega}} [\|\mathbf{q}\|]_h \quad \forall \mathbf{q} \in \mathbf{Q}_h. \qquad (223) $$

In short, by means of the above described boundary interpolation, in the case of mixed boundary conditions we can expect the method's order to become optimal for $k > 0$ too, as long as $\Omega$ is sufficiently smooth, $\mathbf{p} \in [H^{k+1}(\Omega)]^2$, $u \in H^{\max[3,k+1]}(\Omega)$ and $\nabla \cdot \mathbf{p} \in H^{k+1}(\Omega)$. Notice however that a formal assertion in this connection requires the demonstration of the uniform stability of the thus modified formulation. We strongly conjecture that this property holds true, since we can take better advantage of material exploited in Section 5, owing to the enhanced regularity of the underlying auxiliary quantities.

# 8. FINAL REMARKS AND CONCLUSIONS

## 8.1. The case of pure Dirichlet conditions

Although it is not the purpose of our method to handle Dirichlet conditions prescribed on the whole boundary of a curved domain, a comment on this case is in order. More precisely, we mean the problem

$$ \left\{ \begin{array}{ll} \text{Given } f \in L^2(\Omega) \text{ find } (\mathbf{p}; u) \in \mathbf{H}(\mathrm{div}, \Omega) \times L^2(\Omega) \text{ such that} & \\ -(\nabla \cdot \mathbf{p}, v) = (f, v) & \forall v \in L^2(\Omega); \\ (\mathbf{p}, \mathbf{q}) + (u, \nabla \cdot \mathbf{q}) = 0 & \forall \mathbf{q} \in \mathbf{H}(\mathrm{div}, \Omega). \end{array} \right. \qquad (224) $$

Mimicking (7) and observing that we do not need $\mathbf{P}_h$ in this case, (224) can be approximated by

$$ \left\{ \begin{array}{ll} \text{Find } (\mathbf{p}_h; u_h) \in \mathbf{Q}_h \times V_h \text{ such that} & \\ -(\nabla \cdot \mathbf{p}_h, v) = (f, v) & \forall v \in V_h; \\ (\mathbf{p}_h, \mathbf{q}) + (u_h, \nabla \cdot \mathbf{q}) = 0 & \forall \mathbf{q} \in \mathbf{Q}_h. \end{array} \right. \qquad (225) $$

According to the predictions in Section 7, everything works fine for $k = 0$ and $k = 1$, as long as $u \in H^{k+2}(\Omega)$. Otherwise stated, (225) yields optimally converging approximations for $k \le 1$ and furthermore we only have to assume that $\Omega$ is of the $C^k$-class. Indeed, only an inf-sup condition involving the product space $\mathbf{Q}_h \times V_h$ is necessary in this case, i.e., the classical one given in [10]. However, if $k > 1$ the variational residual dominates the interpolation error, and hence the optimal order in the sense of $\mathbf{H}(\mathrm{div}, \Omega_h') \times L^2(\Omega_h')$ can only be attained if we interpolate $\Gamma_T$ by a polynomial arc of degree $k+1$ $\forall T \in S_h$ and integrate in the underlying curved element instead of $T$, as explained in Section 7.3.

## 8.2. Miscellaneous remarks

### 8.2.1. On the assumed geometric regularity of the problem definition domain

The ideal geometric regularity for the error estimates given in Section 7 to hold is the assumption that $\Omega$ be of the piecewise $C^{k+2}$-class (cf. [5]). However, in this work we required much more, though in Lemma 3.3 only. We recall that this result was used as a basis for studying the bilinear form $d_h$ and errors for interpolates in $\mathbf{P}_h$. Actually we believe that this lemma, whose proof heavily relies on properties of Gaussian quadrature in an interval of $\mathbb{R}$, is of academic interest only. Indeed, we obtained results of comparable accuracy with $RT_0$ by prescribing zero normal components on $\Gamma_T$ for all $T \in S_{1,h}$ at points $P$ shifted away from the points located on the perpendicular to $e_T$ through its (Gaussian quadrature) mid-point. As matter of fact, the past experience on the same type of Petrov-Galerkin formulation used in this work, including results reported in [27, 28, 31], indicate that the position of the points on the true boundary at which DOFs are prescribed, does not affect accuracy very much, let alone convergence rates.
### 8.2.2. Petrov–Galerkin approach vs. parametric elements

Some users might wonder if our Pertov–Galerkin approach with different spaces of trial- and test-fields induces higher implementation cost as compared to a pure Galerkin one with parametric elements. Such a conjecture would be due to the need for calculating a different basis for each triangle having two vertexes on $\Gamma_1$. In this respect we should report that two codes to solve Poisson Dirichlet problems in curved domains with both approaches, using Lagrange finite elements based either on the Petrov–Galerkin formulation or on the classical isoparametric formulation were compared in [26]. It turns out that both approaches were equivalent in terms of CPU time, with a slight advantage of the former over the later. On the other hand, the Petrov–Galerkin approach showed up significantly more accurate, which might be due to the fact that exact numerical integration can be used to compute the element matrices, in contrast to the parametric case. Finally the implementation of both approaches requires practically the same amount of data defining the curvilinear boundary. Incidentally, this is certainly true of any other method to handle accurately values prescribed on this kind of manifolds.

### 8.2.3. Possible extension to the three-dimensional case

The natural extension of our method to the three-dimensional case in a context of a similar nature consists of considering the Raviart–Thomas method for tetrahedrons. In principle given normal fluxes have to be prescribed at $dim P_k(F_T) = (k+2)(k+1)/2$ points on $\Gamma_1$ for every face $F_T$ of a tetrahedron in a mesh having three vertexes on this boundary portion. However, in this case some preliminary considerations are in order: A *priori* the definition of the modified Raviart–Thomas space $P_h$ using projections of Gaussian quadrature points onto the curvilinear boundary is not suitable for the three-dimensional case. Indeed, in contrast to the analogous two-dimensional situation exploited here, there is no Gauss quadrature formula with $dim P_k(F_T)$ points yielding exact integrals for polynomials of degree $\le 2k+1$ on a triangular face $F_T$, except for $k=0$. In view of this, a counterpart of Lemma 3.3 must be discarded for arbitrary $k$. Nevertheless, as pointed out in Section 8.2.1, this lemma seems to be only a formal requirement, and DOFs attached to more freely chosen boundary points could be used instead, without downgrading neither accuracy nor convergence rates. But of course such a conjecture must be carefully checked, at first numerically, and eventually confirmed in a formal sense at a subsequent stage. That is why the authors intend to address the three-dimensional case in the near future.

### 8.2.4. Applicability to BDM methods

The same formulation considered here can be applied to BDM methods in two-dimensional space introduced in [11]. In this way the definition of the flux space involving triangles having a curved edge contained on the true boundary could be avoided, although a deeper look into the issue is a must. Notice that the $RT_k$ space strictly contains the space $[P_k]^2$, whose fields can be defined taking only the off-boundary DOFs of $RT_k$. This implies that the way normal fluxes are prescribed on the curved boundary play no role, at least as far as interpolation error estimates in $L^2$ are concerned. In contrast, for BDM triangles the boundary nodes too must be taken into account in order to attain the same optimality. It follows that suboptimal error estimates are expected in this case, unless known normal fluxes are prescribed on the true boundary. In the three-dimensional case (cf. [12]) the problem is even more delicate, for geometric issues related to tetrahedrons having three or even two vertexes on the curved boundary come into play. Such complications explained in [28] rule out working with tetrahedrons having three plane faces and a curved face contained in the boundary. However, we believe that our method can also be advantageously applied to this case, which we intend to examine in due course.

## 8.3. Summary and conclusions

This work is an outstanding application of the theoretical framework built up in [30] to support the reliability analysis of numerical methods based on non symmetric mixed variational formulations with different product spaces for shape- and test-fields. This generalizes the well known theory for symmetric mixed formulations with the same product spaces.

This paper is rather long, since we authors endeavored to study mixed finite element methods for PDEs posed in curved domains, by providing a rigorous treatment of several delicate issues inherent to them, while
keeping in mind accuracy improvement and practical feasibility at a time. This is because, as far as we can see, the simultaneous fulfillment of these three requirements has been overlooked in most works on the subject so far.

As shown in previous articles listed in the bibliography, the type of formulation adopted in this paper lines up with well established alternatives as a valid one to obtain enhanced accuracy in the solution of elliptic equations with prescribed solution values on the boundary of a curved domain. In the case of the mixed formulation of second order equations our method was compared with the formulation employing the same product spaces of shape- and test-fields, defined upon an approximating polygonal domain. For example, in the case of the $RT_1$ method, the mean-square errors of the flux variable are roughly divided by two with the use of our method, while those of the multiplier remain practically the same for both approaches, as much as their observed second order convergence rates in $H(div)$, as reported in [29].

## ACKNOWLEDGMENTS

The authors would like to thank the two referees for their valuable comments and suggestions that greatly contributed to the improvement of this article. They also duly acknowledge the institutional support for enabling the successful accomplishment of their work.

## DATA AVAILABILITY STATEMENT

The research data associated with this article are included in the article.

## SUPPLEMENTARY MATERIAL

The Supplementary Material is available at https://www.esaim-m2an.org/10.1051/m2an/2025028/olm.

## REFERENCES

[1] R.A. Adams, Sobolev Spaces. Academic Press, New York (1975).

[2] I. Babuška, Error bound for the finite element method. *Num. Math.* **16** (1971) 322–333.

[3] I. Babuška, The finite element method with Lagrangian multipliers. *Num. Math.* **20** (1973) 179–192.

[4] F. Bertrand and V. Ruas, A variant of the Raviart–Thomas method for smooth domains using straight-edged triangles. Preprint arXiv:2307.03503 (2025).

[5] F. Bertrand and G. Starke, Parametric Raviart–Thomas mixed methods on domains with curved surfaces. *SIAM J. Numer. Anal.* **54** (2016) 3648–3667.

[6] F. Bertrand, S. Münzenmaier and G. Starke, First-order system least-squares on curved boundaries: higher-order Raviart–Thomas elements. *SIAM J. Num. Anal.* **52** (2014) 3165–3180.

[7] E.K. Blum, Numerical Analysis and Computation: Theory and Practice. Addison Wesley (1972).

[8] S.C. Brenner and L.R. Scott, The Mathematical Theory of Finite Element Methods. *Texts Appl. Math.* Vol. 15. Springer (2008).

[9] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers. *RAIRO Anal. Numér.* **8** (1974) 129–151.

[10] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods. Springer-Verlag (1991).

[11] F. Brezzi, J. Douglas Jr. and D. Marini, Two families of mixed finite elements for second order elliptic problems. *Numer. Math.* **47** (1985) 217–235.

[12] F. Brezzi, J. Douglas Jr., R. Durán and M. Fortin, Mixed finite elements for second order elliptic problems in three variables. *Numer. Math.* **51** (1987) 237–250.

[13] P.G. Ciarlet, The Finite Element Method for Elliptic Problems. North Holland (1978).

[14] P.G. Ciarlet, Three-Dimensional Elasticity. Vol. 1. Elsevier Science (1994).

[15] P.G. Ciarlet and P.A. Raviart, The combined effect of curved boundaries and numerical integration in isoparametric finite element methods, in *The Mathematical Foundations of the Finite Element Method with Applications to Part. Diff. Eqns.*, edited by A.K. Aziz. Academic Press (1972) 409–474.

[16] J.A. Cuminato and V. Ruas, Unification of distance inequalities for linear variational problems. *Comput. Appl. Math.* **34** (2015) 1009–1033.
[17] E. Di Nezza, G. Palatucci and E. Valdinoci, Hitchhiker's guide to the fractional Sobolev spaces. *Bull. Sci. Math.* **136** (2012) 521–573.

[18] G. Duvaut and J.-L. Lions, Les Inéquations en Mécanique et en Physique. Dunod (1972).

[19] A. Ern and J.L. Guermond, Finite Elements I, Approximation and Interpolation. Springer (2021).

[20] E.H. Georgoulis, Inverse-type estimates on hp-finite element spaces and applications. *Math. Comput.* **77** (2008) 201–219.

[21] V. Girault and P.A. Raviart, Finite Element Methods for Navier-Stokes Equations: Theory and Algorithms. Springer (1986).

[22] P. Grisvard, Elliptic Problems in Non Smooth Domains. Pitman Publ. Inc. (1985).

[23] A.J. Lew and M. Negri, Optimal convergence of a discontinuous-Galerkin-based immersed boundary method. *ESAIM M2AN* **45** (2011) 251–274.

[24] A. Quarteroni, R. Sacco and F. Saleri, Numerical Mathematics. Texts Appl. Math. Springer (2007).

[25] P.-A. Raviart and J.-M. Thomas, Mixed finite element methods for second order elliptic problems, in Lecture Notes in Mathematics. Vol. 606. Springer Verlag (1977) 292–315.

[26] V. Ruas, Optimal simplex finite-element approximations of arbitrary order in curved domains circumventing the isoparametric technique. Preprint arXiv:1701.00663 (2017).

[27] V. Ruas, Optimal Lagrange and Hermite finite elements for Dirichlet problems in curved domains with straight-edged triangles. *ZAMM* **100** (2020) e-2019002962020.

[28] V. Ruas, Optimal-rate nonconforming finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra. *IMA J. Num. Anal.* **41** (2021) 1368–1410.

[29] V. Ruas, Success and failure of attempts to improve the accuracy of Raviart-Thomas mixed finite elements in curved domains. *Examples Counterexamples* **2** (2022) 100037.

[30] V. Ruas, Fine error bounds for approximate asymmetric saddle point problems. *Comput. Appl. Math.* **43** (2024) 160.

[31] V. Ruas and M.A. Silva Ramos, Efficiency on non parametric finite elements for optimal-order enforcement of Dirichlet conditions on curvilinear boundaries. *J. Comput. Appl. Math.* **394** (2021) 113523.

[32] A.M. Sanchez and R. Arcangeli, Estimations des erreurs de meilleure approximation polynomiale et d'interpolation de Lagrange dans les espaces de Sobolev d'ordre non entier. *Numer. Math.* **45** (1984) 301–321.

[33] D.B. Stein, R.D. Guy and B. Thomases, Immersed boundary smooth extension: a high-order method for solving PDE on arbitrary smooth domains using Fourier spectral methods. *J. Comput. Phys.* **304** (2016) 252–274.

[34] G. Strang and G. Fix, An Analysis of the Finite Element Method. Prentice-Hall (1973).

[35] J.L. Taylor, S. Kim and R.M. Brown, The Green function for elliptic systems in two dimensions. *Comm. Part. Differ. Eqns.* **38** (2013) 1574–1600.

[36] R. Verfürth, A Posteriori Error Estimation Techniques for Finite Element Methods. Oxford Sci. Publ. (2013).

[37] O.C. Zienkiewicz, The Finite Element Method in Engineering Science. McGraw-Hill (1971).

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

Please help to maintain this journal in open access!

This journal is currently published in open access under the **Subscribe to Open model (S2O)**. We are thankful to our subscribers and supporters for making it possible to publish this journal in open access in the current year, free of charge for authors and readers.

Check with your library that it subscribes to the journal, or consider making a personal donation to the S2O programme by contacting subscribers@edpsciences.org.

More information, including a list of supporters and financial transparency reports, is available at https://edpsciences.org/en/subscribe-to-open-s2o.