m}|_{k+1,\tilde{T}}^{2}+|\nabla\cdot\tilde{\bf m}|_{k+1,T^{\Gamma}}^{2}\\}.$
(203)

Assuming that $h_{T}\leq 1/(8C_{I3})^{2}$, we combine (203) and (184) to
conclude that the interpolation error of $\nabla\cdot\tilde{\bf m}$ can be
estimated by

$\parallel\nabla\cdot[\tilde{\bf m}-\tilde{\Pi}_{h}\tilde{\bf
m}]\parallel_{0,h}\leq C_{k,div}h^{k+1}\\{\|\tilde{\bf
m}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}+\|\tilde{\bf
m}\|_{k+1,\tilde{\Omega}}+|\nabla\cdot\tilde{\bf m}|_{k+1,\tilde{\Omega}}\\},$
(204)

for a suitable constant $C_{k,div}$ independent of $h$.
Owing to (204) and (183), the result immediately follows with
$C_{H}=[C_{k,div}^{2}+C_{k,\tilde{\Omega}}^{2}]^{1/2}$

##  6 Error estimates

In this section we endeavor to obtain error estimates for problems (10) and
(115).
In the case of curvilinear boundaries, in principle, we are dealing with an
external approximation also called a variational crime (cf. [30]). Hence we
resort to the general upper bound for linear variational problems given in
[13] for weakly coercive bilinear forms defined upon two different spaces.

###  6.1 The convex case

In order to simplify things a little we first address the case where $\Omega$
is convex.
To begin with, we extend the bilinear form $c:({\bf P}_{h,k}\times
V_{h,k})\times({\bf Q}_{h,k}\times V_{h,k})\rightarrow\Re$ defined by (23) to
$({\bf Q}\times V)\times({\bf Q}\times V)$. Since we are assuming that
$\Omega$ is convex, every pair $({\bf r};w)\in{\bf Q}\times V$ is defined in
the whole $\Omega_{h}$ and therefore such an extension of $c$ can be expressed
exactly as in (23). We recall that the bilinear form $c$ is bounded in $({\bf
Q}\times V)\times({\bf Q}\times V)$ with norm $\|c\|$ (cf. Section 4), that
is,

$\ c(({\bf r};w),({\bf q};v))\leq\|c\||||({\bf r};w)|||_{h}|||{\bf
q};v)|||_{h}\;\forall({\bf r};w)\in({\bf Q}+{\bf P}_{h,k})\times
V,\;\forall({\bf q};v)\in{\bf Q}_{h,k}\times V_{h,k}.$ (205)

We next consider separately two possibilities, namely,

I) $\Omega$ is convex and $length(\Gamma_{0})>0$
In view of the theoretical framework supplied in [13], recalling (113), as
long as $h$ is sufficiently small and $\Omega$ is smooth enough, we can write

$\left\\{\begin{array}[]{l}|||({\bf p};u)-({\bf
p}_{h};u_{h})|||_{h}\leq\displaystyle\frac{1}{\gamma}\left[\|c\|\inf_{({\bf
r};w)\in{\bf P}_{h,k}\times V_{h,k}}|||({\bf p};u)-({\bf
r};w)|||_{h}\right.\\\ \left.+\displaystyle\sup_{({\bf q};v)\in{\bf
Q}_{h,k}\times V_{h,k}/|||({\bf q};v)|||_{h}=1}|c(({\bf p};u),({\bf
q};v))-(f,v)_{h}|.\right]\end{array}\right.$ (206)

The $inf$ term in (206) can be bounded above by $|||({\bf
p};u)-(\tilde{\Pi}_{h}{\bf p};\tilde{\pi}_{h}u)|||_{h}$, where
$\tilde{\pi}_{h}u$ is the standard Lagrange interpolate of $u$ for $k\geq 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,\ldots$ (cf. [11])

$\|u-\tilde{\pi}_{h}u\|_{0,h}\leq C_{\pi}h^{k+1}|u|_{k+1}.$ (207)

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

${\mathcal{R}}_{h,k}(u,{\bf q}):=({\bf p},{\bf q})_{h}-(u,\nabla\cdot{\bf
q})_{h}=-(u,{\bf q}\cdot{\bf n}_{h})_{0,\Gamma_{0,h}}\;\forall{\bf q}\in{\bf
Q}_{h,k},$ (208)

In principle, ${\mathcal{R}}_{h,k}(u,{\bf 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 6.1

Let $length(\Gamma_{0})>0$ and ${\mathcal{R}}_{h,k}(u,{\bf q})$ be given by
(208). If the first (resp. second) component ${\bf p}$ (resp. u) of the
solution of (2) 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\geq 0$ it holds

$|{\mathcal{R}}_{h,k}(u,{\bf q})|\leq C_{\Gamma_{0}}h^{3/2}\|u\|_{2}[|{\bf
q}|]_{h}\;\forall{\bf q}\in{\bf Q}_{h,k}.$ (209)

Proof. First of all we have

$(u,{\bf q}\cdot{\bf n}_{h})_{e_{T}}\leq\|u\|_{0,e_{T}}\|{\bf q}\cdot{\bf
n}_{h}\|_{0,e_{T}}\;\forall T\in{\mathcal{S}}_{0,h}.$ (210)

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

$\|{\bf q}\cdot{\bf n}_{h}\|_{0,e_{T}}\leq l_{T}^{1/2}\|{\bf
q}\|_{0,\infty,T}\leq C_{0I}h_{T}^{-1/2}\|{\bf q}\|_{0,T}.$ (211)

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}}\leq C_{e0}h_{T}^{j}\|\partial^{j}u/\partial
x^{j}\|_{0,e_{T}}\mbox{ for }1\leq j\leq 2.$ (212)

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.$ (213)

${\bf t}$ being the tangent unit vector along $\Gamma$ we introduce the direct
reference frame ${\mathcal{R}}:=(N;{\bf t},{\bf n})$. Then recalling the angle
$\varphi$ between ${\bf n}_{h}$ and ${\bf n}$, let the unit vector along
$e_{T}$ be ${\bf l}_{T}$. We have ${\bf l}_{T}=\cos{\varphi}{\bf
t}-\sin{\varphi}{\bf n}$ and hence we have $\partial u/\partial x=\nabla
u_{|{\mathcal{R}}}\cdot(\cos{\varphi};-\sin{\varphi})$. We know that
$|\sin{\varphi}|\leq C_{\Gamma}h_{T}$ owing to (3) 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))|\leq C_{\Gamma}h_{T}|(\partial u/\partial
n)(x,\phi(x))|\;\forall x\in[0,l_{T}].$ (214)

Plugging (214) into (213) we obtain

$|(\partial u/\partial x)(x,0)|\leq C_{\Gamma}h_{T}|(\partial u/\partial
n)(x,\phi(x))|+C_{0}^{1/2}h_{T}\displaystyle\left[\int_{M}^{N}|\partial^{2}u/\partial
x\partial y|^{2}dy\right]^{1/2}\;\forall x\in[0,l_{T}].$ (215)

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 (215) into (212), it trivially holds for a certain
constant $C_{e1}$ independent of $T$:

$\|u\|^{2}_{0,e_{T}}\leq C^{2}_{e1}h_{T}^{4}[\|\partial u/\partial
n\|^{2}_{0,\Gamma_{T}}+|u|^{2}_{2,\Delta_{T}}]$ (216)

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

$(u,{\bf q}\cdot{\bf n}_{h})_{0,\Gamma_{0,h}}\leq
C_{0I}C_{e1}h^{3/2}\\{\|\partial u/\partial
n\|_{0,\Gamma_{0}}+|u|_{2}\\}[|{\bf q}|]_{h}.$ (217)

Finally, thanks to the Trace Theorem for $H^{2}(\Omega)$ (cf. [1]), (209)
derives from (217).

As a consequence of Proposition 6.1 we trivially have

###### Theorem 6.2

Assume that $length(\Gamma_{0})>0$, $\Omega$ is smooth enough, $h$ is
sufficiently small and $k\geq 0$. If for a small $\varepsilon>0$ $\tilde{\bf
p}\in[H^{k+1+\varepsilon(1-\min[1,k])}(\tilde{\Omega})]^{2}$, $\nabla\cdot{\bf
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) by the solution of (10).

$|||({\bf p};u)-({\bf p}_{h};u_{h})|||_{h}\leq\displaystyle
C_{01}\left\\{h^{3/2}\|u\|_{2}+h^{k+1}\left[\|\tilde{\bf
p}\|_{k+1+\varepsilon(1-\min[1,k]),\tilde{\Omega}}+|\nabla\cdot{\bf
p}|_{k+1}+|u|_{k+1}\right]\right\\}.\mbox{$\rule{5.69054pt}{5.69054pt}$}$
(218)

###### Remark 1

According to the estimate (218), 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)$ [19]. Therefore
${\bf 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 Subsection 6.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$
Since $u\in L^{2}_{0}(\Omega)$ in this case, a natural way to define the
interpolate of $u$ in $V_{h,k}$ is as follows: We start with the operator
$\tilde{\pi}_{h}$ from $V^{k}(\Omega):=H^{k+1}(\Omega)\cap V$ into
$L^{2}(\Omega_{h})$ already defined in I). Noticing that $\tilde{\pi}_{h}u$
does not necessarily belong to $L^{2}_{0}(\Omega_{h})$, $\tilde{\pi}_{h}u$ may
not belong to $V_{h,k}$. Thus we set
$\bar{\pi}_{h}u:=\tilde{\pi}_{h}u+C_{h}\in V_{h,k}$ with
$C_{h}:=-\int_{\Omega_{h}}\tilde{\pi}_{h}u/area(\Omega_{h})$. While on the one
hand the estimate of the interpolation error in $V_{h,k}$ of $u$ is not
optimal, we note that the estimate $\|u-\tilde{\pi}_{h}u\|_{0,h}\leq
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 ${\bf P}_{h,k}\times V_{h,k}$ of the pair $({\bf p};\bar{u})$
instead of $({\bf p};u)$, as seen below.
It is noticeable that, since $\int_{\Omega}\nabla\cdot{\bf q}=0$ $\forall{\bf
q}\in{\bf Q}$, the replacement of $u$ by $\bar{u}$ does not change the essence
of problem (2); it is just a way like any other to fix the additive constant
up to which $u$ is defined.
Summarizing, the following error bound applies to the pair $({\bf
p},\bar{u})$.

$\left\\{\begin{array}[]{l}|||({\bf p};\bar{u})-({\bf
p}_{h};u_{h})|||_{h}\leq\displaystyle\frac{1}{\gamma}\left[\|c\|\inf_{({\bf
r};w)\in{\bf P}_{h,k}\times V_{h,k}}|||({\bf p};\bar{u})-({\bf
r};w)|||_{h}\right.\\\ \left.+\displaystyle\sup_{({\bf q};v)\in{\bf
Q}_{h,k}\times V_{h,k}/|||({\bf q};v)|||_{h}=1}|c(({\bf p};\bar{u}),({\bf
q};v))-(f,v)_{h}|.\right]\end{array}\right.$ (219)

It turns out that

$\displaystyle\inf_{w\in
V_{h,k}}\|\bar{u}-w\|_{0,h}\leq\|\bar{u}-\bar{\pi}_{h}u\|_{0,h}=\|u-\tilde{\pi}_{h}u\|_{0,h}\leq
C_{\pi}h^{k+1}|u|_{k+1}.$ (220)

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

$\inf_{({\bf r};w)\in{\bf P}_{h,k}\times V_{h,k}}|||({\bf p};\bar{u})-({\bf
r};w)|||_{h}\leq C_{I}h^{k+1}\displaystyle\left\\{\|\tilde{\bf
p}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}+|\nabla\cdot{\bf
p}|_{k+1}+|u|_{k+1}\right\\}.$ (221)

As for the $sup$ term in (206), we can assert that it vanishes identically if
$length(\Gamma_{0})=0$. Indeed, $\nabla\cdot{\bf p}=f$ in $\Omega_{h}$ and on
the other hand, as long as $u\in H^{k+1}(\Omega)$,

$({\bf p},{\bf q})_{h}-(\bar{u},\nabla\cdot{\bf q})_{h}=({\bf
p}-\nabla\bar{u},{\bf q})_{h}=0\;\forall{\bf q}\in{\bf Q}_{h,k},$ (222)

since ${\bf q}\cdot{\bf n}_{h}$ vanishes everywhere in $\Gamma_{h}$ in this
case.
This leads to the following result.

###### Theorem 6.3

Assume that $length(\Gamma_{0})=0$, $\Omega$ is smooth enough and $h$ is
sufficiently small. If $\tilde{\bf
p}\in[H^{k+1+\varepsilon(1-\min[k,1])}(\tilde{\Omega})]^{2}$, $\nabla\cdot{\bf
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) by the solution of
(10).

$\left\\{\begin{array}[]{l}|||({\bf p};\bar{u})-({\bf
p}_{h};u_{h})|||_{h}\leq\displaystyle\frac{C_{I}\|c\|}{\gamma}h^{k+1}\left\\{\|\tilde{\bf
p}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}+|\nabla\cdot{\bf
p}|_{k+1}+|u|_{k+1}\right\\}\\\ \mbox{where }u-\bar{u}=\mbox{ constant, for
}k=0,1,2,\ldots.\mbox{$\rule{5.69054pt}{5.69054pt}$}\end{array}\right.$ (223)

###  6.2 The non convex case

New issues must be overcome in the non convex case, as far as error estimates
are concerned.
First of all, if $\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\bar{\Omega}\neq\emptyset$. In order to remedy such a
shortcoming, akin to [24], we use rather the norms related to the domain
$\Omega_{h}^{{}^{\prime}}$ instead of $\Omega_{h}$, in accordance with the
definitions and notations given in Section 2. We also refer to Subsection 4.2
for some definitions and concepts evoked hereafter. Here again, we work with a
smooth open set $\tilde{\Omega}$ sufficiently large to strictly contain not
only $\Omega$, but also the sets $\tilde{\Omega}_{h}$ for all
${\mathcal{T}}_{h}\in{\mathcal{F}}^{{}^{\prime}}$.
Further, for all meshes of the family in use, it is convenient to (abstractly)
extend somehow outside $\Omega$ the exact solution $({\bf p};u)$ of (1). Let
$(\tilde{\bf p};\tilde{u})$ be such an extension of $({\bf p};u)$ defined in
$\tilde{\Omega}$. The pair $(\tilde{\bf p};\tilde{u})$ must have the required
regularity for method’s order to remain the same as in the convex case. Akin
to [24], the construction of the extension $\tilde{u}$ of $u$ \- which
incidentally induces a regular extension $\tilde{f}$ of $f$ in the form of
$-\Delta\tilde{u}$ -, is the one advocated in [29]. Naturally enough we set
$\tilde{\bf p}:=\nabla\tilde{u}$. Of course, the explicit knowledge of neither
$(\tilde{\bf p};\tilde{u})$ nor $\tilde{f}$ is necessary. The point here is
that all these quantities do exist, and can thus be employed in support to the
theoretical reliability analysis of the non convex case, which we next
complete.

First notice that the extension technique allows for establishing an estimate
for the interpolation error in the non convex case analogous to (182). More
precisely if $\tilde{\bf p}\in{\mathcal{M}}^{k}(\tilde{\Omega})$ we have

$[|\tilde{\bf p}-\tilde{\Pi}_{h}\tilde{\bf p}|]_{h}\leq
C_{H}h^{k+1}[\|\tilde{\bf
p}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}^{2}+|\nabla\cdot\tilde{\bf
p}|_{k+1,\tilde{\Omega}}^{2}]^{1/2}.$ (224)

Furthermore, we may exploit the extension $\tilde{f}$ in the framework of the
approximate right hand side $F_{h}(v)$ of problem (115). 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\Omega$ $\forall
T\in{\mathcal{S}}_{h}$ for all ${\mathcal{T}}_{h}$ in
${\mathcal{F}}^{{}^{\prime}}$.
Then we set

$F_{h}(v):=(f_{h},v)_{h}\;\forall v\in V_{h,k}.$ (225)

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

$\|\tilde{f}-f_{h}\|_{0,h}\leq
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 (225)

$|(f,v)_{h}-F_{h}(v)|\leq
C_{\chi,k}h^{k+1}\|\tilde{f}\|_{k+1+\varepsilon(1-\min[k,1]),\Omega_{h}}\|v\|_{0,h}\forall
v\in V_{h,k}$ (226)

Henceforth we definitively take $\tilde{f}=-\Delta\tilde{u}$ in
$\tilde{\Omega}\setminus\Omega$ as the extension of $f$ with the required
smoothness.
Notice that, thanks to (113) with $\gamma=\gamma^{{}^{\prime}}$ given by
(178), we can assert that the approximate problem (115) has a unique solution,
and moreover that the constant $\gamma^{{}^{\prime}}>0$ is independent of $h$.
It is noteworthy that, in the non convex case, the bilinear form $c$ has to be
extended to pairs in ${\bf H}(div,\Omega_{h})\times L^{2}(\Omega_{h})$, which
is not the case of $({\bf p};u)$, but rather of $(\tilde{\bf p};\tilde{u})$ by
construction. Incidentally, this extension does not modify at all the value of
$\|c\|$.

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

III) $\Omega$ is not convex and $length(\Gamma_{0})>0$
First of all, in perfect analogy to (206), the following upper bound applies
to the approximation error:

$\left\\{\begin{array}[]{l}|||({\bf p};u)-({\bf
p}_{h};u_{h})|||_{h}^{{}^{\prime}}\leq|||(\tilde{\bf p};\tilde{u})-({\bf
p}_{h};u_{h})|||_{h}\leq\displaystyle\frac{1}{\gamma^{{}^{\prime}}}\left[\|c\|\inf_{({\bf
r};w)\in{\bf P}_{h,k}\times V_{h,k}}|||(\tilde{\bf p};\tilde{u})-({\bf
r};w)|||_{h}\right.\\\ \left.+\displaystyle\sup_{({\bf q};v)\in{\bf
Q}_{h,k}\times V_{h,k}/|||({\bf q};v)|||_{h}=1}|c((\tilde{\bf
p};\tilde{u}),({\bf q};v))-F_{h}(v)|.\right]\end{array}\right.$ (227)

We refrain from giving a detailed justification of inequality (227) since it
derives from well known arguments (see e.g. [24]).
We next focus on the estimation of the contribution to the approximation error
of the variational residual incorporating the functional $F_{h}(v)$ instead of
$(\tilde{f},v)_{h}$, namely, $|c((\tilde{\bf p};\tilde{u}),({\bf
q};v))-F_{h}(v)|$. With this aim we first write

$|c((\tilde{\bf p};\tilde{u}),({\bf q};v))-F_{h}(v)|\leq|c((\tilde{\bf
p};\tilde{u}),({\bf q};v))-(\tilde{f},v)_{h}|+|(\tilde{f},v)_{h}-F_{h}(v)|.$
(228)

Now we observe that, either by construction, or due to (1), it holds
$-\nabla\cdot\tilde{\bf p}=\tilde{f}$ and $\tilde{\bf 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 (228) is
${\mathcal{R}}_{h,k}(\tilde{u},{\bf q}):=(\tilde{u},{\bf q}\cdot{\bf
n}_{h})_{0,\Gamma_{0,h}}$. It is easy to see that
${\mathcal{R}}_{h,k}(\tilde{u},{\bf q})$ can be estimated in the same manner
as in (209), that is,

$|{\mathcal{R}}_{h,k}(\tilde{u},{\bf q})|\leq
C_{\Gamma_{0}}h^{3/2}\|\tilde{u}\|_{2,\tilde{\Omega}_{h}}[|{\bf
q}|]_{h}\;\forall{\bf q}\in{\bf Q}_{h,k}.$ (229)

Combining (226) and (229) and setting
$\tilde{C}_{r}:=\max[C_{\Gamma_{0}},C_{\chi,k}]$, we immediately come up with

$\displaystyle\sup_{({\bf q};v)\in{\bf Q}_{h,k}\times V_{h,k}/|||({\bf
q};v)|||_{h}=1}|c((\tilde{\bf p};\tilde{u}),({\bf
q};v))-F_{h}(v)|\leq\tilde{C}_{r}\displaystyle\left[h^{\frac{3}{2}}\\!\|\tilde{u}\|_{2,\tilde{\Omega}}+h^{k+1}\|\tilde{f}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}\right].$
(230)

Now taking $\tilde{\bf m}=\tilde{\bf p}$ in (224) and letting
$\tilde{\pi}_{h}\tilde{u}$ be the $V_{h,k}$-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:

$\begin{array}[]{l}\inf_{({\bf r};w)\in{\bf P}_{h,k}\times
V_{h,k}}|||(\tilde{\bf p};\tilde{u})-({\bf r};w)|||_{h}\leq|||(\tilde{\bf
p};\tilde{u})-(\tilde{\Pi}_{h}\tilde{\bf
p};\tilde{\pi}_{h}\tilde{u})|||_{h}\\\
\leq\tilde{C}_{I}h^{k+1}\displaystyle\left\\{\|\tilde{\bf
p}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}+|\nabla\cdot\tilde{\bf
p}|_{k+1,\tilde{\Omega}}+|\tilde{u}|_{k+1,\tilde{\Omega}}\right\\}.\end{array}$
(231)

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

###### Theorem 6.4

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

$\begin{array}[]{l}|||({\bf p};u)-({\bf
p}_{h};u_{h})|||_{h}^{{}^{\prime}}\leq\max[\tilde{C}_{I},\tilde{C}_{r}]\displaystyle\left\\{h^{3/2}\|\tilde{u}\|_{2,\tilde{\Omega}}\right.\\\
\left.+h^{k+1}\displaystyle\left[\|\tilde{f}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}+\|\tilde{\bf
p}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}+|\nabla\cdot\tilde{\bf
p}|_{k+1,\tilde{\Omega}}+|\tilde{u}|_{k+1,\tilde{\Omega}}\right]\right\\}.\mbox{
$\rule{5.69054pt}{5.69054pt}$}\end{array}$ (232)

From (232) 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,k}(\tilde{u},{\bf
q})|\leq\|\tilde{u}\|_{1/2,\Gamma_{0,h}}\|{\bf q}\cdot{\bf
n}_{h}\|_{-1/2,\Gamma_{0,h}},$ (233)

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 [18] there is a
constant $C_{\Gamma_{0},h}$ depending on $\Omega_{h}$ such that

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

On the other hand we have

$\|\tilde{u}\|_{1/2,\Gamma_{0,h}}=\displaystyle\left[\sum_{T\in{\mathcal{S}}_{0,h}}\|\tilde{u}\|_{1/2,e_{T}}^{2}\right]^{1/2}.$
(235)

It turns out that $\tilde{u}$ vanishes at the end points of $e_{T}$.
Therefore, by standard approximation results (cf. [28]) 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}}.$ (236)

Plugging (236) into (235) and combining the resulting expression with (234),
after strightforward calculations we conclude from (233) that

$|{\mathcal{R}}_{h,k}(\tilde{u},{\bf q})|\leq
C_{u}C_{\Gamma_{0},h}h^{2}|\tilde{u}|_{5/2,\Gamma_{0,h}}[|{\bf q}|]_{h}.$
(237)

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}}.$ (238)

Plugging (238) into (237) we conclude that the residual
${\mathcal{R}}_{h,k}(\tilde{u},{\bf 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,k}(\tilde{u},{\bf q})|\leq
C_{\tilde{\Omega}_{h}}h^{2}\|\tilde{u}\|_{3,\tilde{\Omega}}[|{\bf q}|]_{h}.$
(239)

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 the next section we will 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$
To begin with we note that a modification analogous to (219) applies. 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}/area(\Omega_{h})$. Then we
note that the estimate $\|\tilde{u}-\tilde{\pi}_{h}\tilde{u}\|_{0,h}\leq
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,k}$ and thus

$\left\\{\begin{array}[]{l}\|\bar{u}-\bar{\pi}_{h}\tilde{u}\|_{0,h}\leq
C_{\pi}h^{k+1}|\tilde{u}|_{k+1,\tilde{\Omega}}\\\ \mbox{where}\\\
\bar{u}=\tilde{u}+C_{h}\mbox{ and
}\bar{\pi}_{h}\tilde{u}:=\tilde{\pi}_{h}\tilde{u}+C_{h}.\\\ \mbox{with
}C_{h}:=-\int_{\Omega_{h}}\tilde{\pi}_{h}\tilde{u}/area(\Omega_{h}).\end{array}\right.$
(240)

(240) can be used to estimate the term

$\displaystyle\inf_{w\in
V_{h,k}}\|\bar{u}-w\|_{0,h}\leq\|\bar{u}-\bar{\pi}_{h}\tilde{u}\|_{0,h}=\|\tilde{u}-\tilde{\pi}_{h}\tilde{u}\|_{0,h}\leq
C_{\pi}h^{k+1}|\tilde{u}|_{k+1,\tilde{\Omega}}$ (241)

in the following error bound that holds for problem (115):

$\left\\{\begin{array}[]{l}|||({\bf p};\bar{u})-({\bf
p}_{h};u_{h})|||^{{}^{\prime}}_{h}\leq\displaystyle\frac{1}{\gamma^{{}^{\prime}}}\left[\|c\|\inf_{({\bf
r};w)\in{\bf P}_{h,k}\times V_{h,k}}|||(\tilde{\bf p};\bar{u})-({\bf
r};w)|||_{h}\right.\\\ \left.+\displaystyle\sup_{({\bf q};v)\in{\bf
Q}_{h,k}\times V_{h,k}/|||({\bf q};v)|||_{h}=1}|c((\tilde{\bf
p};\tilde{u}),({\bf q};v))-F_{h}(v)|.\right]\\\ \mbox{where
}\bar{u}=\tilde{u}+C_{h}.\end{array}\right.$ (242)

Rigorously we should have written $\bar{u}$ instead of $\tilde{u}$ in the sup-
term of (242). However, it is legitimate to replace $\bar{u}$ by $\tilde{u}$
therein since $(C_{h},\nabla\cdot{\bf q})_{h}=0$ $\forall{\bf q}\in{\bf
Q}_{h,k}$ if $length(\Gamma_{0})=0$.
Recalling (219), in case $length(\Gamma_{0})=0$ an error estimate like (232)
holds without the term in $h^{3/2}$, though for a pair $({\bf p};\bar{u})$
instead of $({\bf 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 (242). Thus this
leads to

###### Theorem 6.5

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

$\begin{array}[]{l}|||({\bf p};\bar{u})-({\bf
p}_{h};u_{h})|||_{h}^{{}^{\prime}}\leq\max[\tilde{C}_{I},C_{\chi,h}]h^{k+1}\\\
\times\displaystyle\left[\|\tilde{f}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}+\|\tilde{\bf
p}\|_{k+1+\varepsilon(1-\min[k,1]),\tilde{\Omega}}+|\nabla\cdot\tilde{\bf
p}|_{k+1,\tilde{\Omega}}+|\tilde{u}|_{k+1,\tilde{\Omega}}\right].\mbox{$\rule{5.69054pt}{5.69054pt}$}\end{array}$
(243)

Theorem 6.5 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\geq 0$, provided the exact solution of (2) is
sufficiently smooth.

###### Remark 2

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
$\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}^{{}^{\prime}}$ \-
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.

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

As pointed out in Remark 1 of Subsection 6.1, the error estimate (218)
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 Subsection 6.1, the problem stems from the term in the
variational residual involving an integral along $\Gamma_{0,h}$, which is an
${\mathcal{O}}(h^{3/2})$ whatever $k$. Nevertheless this issue can be settled
for $k>0$ if every (straight) triangle in $T\in{\mathcal{S}}_{0,h}$ is
replaced by a triangle $\breve{T}$ having the same straight edges as $T$
except $e_{T}$, which is replaced by a curvilinear edge $\breve{e}_{T}$ having
the same end-points as $e_{T}$. The idea here is not to use a parametric
representation, since the master element $\widehat{\breve{T}}$ associated with
$\breve{T}$ is also a triangle with a curvilinear edge. Actually
$\widehat{\breve{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 $\breve{T}$ itself or in $\widehat{\breve{T}}$.
The point here is that the curvilinear edge $\breve{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{\bf q}\cdot{\bf n}_{h}$ on $\Gamma_{0,h}$, with
${\bf q}\in{\bf Q}_{h,k}$, becomes

$\breve{\mathcal{R}}_{h,k}(\tilde{u},{\bf
q}):=\displaystyle\sum_{T\in{\mathcal{S}}_{0,h}}(\tilde{u},{\bf
q}\cdot\breve{\bf n}_{T})_{\breve{e}_{T}},$ (244)

where ${\bf q}$ is eventually extended to $\breve{T}\setminus T$ in case this
set is not empty and $\breve{\bf n}_{T}$ is the outer normal vector to
$\breve{e}_{T}$.
Likewise the convex case (for $k\geq 1$), we assume that $u\in H^{3}(\Omega)$,
or yet that $\tilde{u}\in H^{3}(\tilde{\Omega})$.
Owing to the construction of $\breve{e}_{T}$,
$\breve{\mathcal{R}}_{h,k}(\tilde{u},{\bf q})$ can be estimated as an
${\mathcal{O}}(h^{k+3/2})$, as long as $\tilde{u}\in H^{3}(\tilde{\Omega})$.
Let us give a precise justification of such an assertion.
First of all, similarly to (210)-(211), we can establish the existence of a
constant $\breve{C}_{0I}$ independent of $T$ such that

$(\tilde{u},{\bf q}\cdot\breve{\bf
n}_{h})_{\breve{e}_{T}}\leq\breve{C}_{0I}h_{T}^{-1/2}\|\tilde{u}\|_{0,\breve{e}_{T}}\|{\bf
q}\|_{0,T}\;\forall T\in{\mathcal{S}}_{0,h}.$ (245)

Therefore, recalling (244) we have

$\breve{\mathcal{R}}_{h,k}(\tilde{u},{\bf
q})\leq\breve{C}_{0I}\displaystyle\left[\sum_{T\in{\mathcal{S}}_{0,h}}h_{T}^{-1}\|\tilde{u}\|^{2}_{0,\breve{e}_{T}}\right]^{1/2}[|{\bf
q}|]_{h}.$ (246)

Now let $P$ be the point of $\breve{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)=\displaystyle\int_{N}^{P}\frac{\partial\tilde{u}}{\partial
y}\;dy\;\forall P\in\breve{e}_{T}.$ (247)

Let $\phi_{h}(x)$ be the function that uniquely represents $\breve{e}_{T}$ in
such a frame. The curvilinear abscissa $s_{h}$ along $\breve{e}_{T}$ satisfies
$ds_{h}=\sqrt{1+[\phi_{h}^{{}^{\prime}}(x)]^{2}}dx$. On the other hand, by
basic interpolation theory in one-dimensional space (see e.g. [6], [21]),
there exists a constant $C_{L}$ independent of $l_{T}$ such that
$\|\phi_{h}^{{}^{\prime}}\|_{0,\infty,e_{T}}\leq
C_{L}\|\phi^{{}^{\prime}}\|_{0,\infty,e_{T}}\leq C_{\Gamma}C_{L}h_{T}$. Hence
from (247) we easily establish the existence of a constant $\breve{C}_{L}$
independent of $T$ fulfilling

$\int_{\breve{e}_{T}}\tilde{u}^{2}(P)\;ds_{h}=\displaystyle\int_{\breve{e}_{T}}\left[\int_{N}^{P}\frac{\partial\tilde{u}}{\partial
y}\;dy\right]^{2}ds_{h}\leq\breve{C}_{L}\displaystyle\int_{0}^{l_{T}}\left[\int_{N}^{P}\frac{\partial\tilde{u}}{\partial
y}\;dy\right]^{2}dx.$ (248)

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

$h_{T}^{-1}\|\tilde{u}\|_{0,\breve{e}_{T}}^{2}\leq\breve{C}_{L}\max_{x\in[0,l_{T}]}|\phi_{h}(x)-\phi(x)|^{2}\|\tilde{u}\|_{1,\infty,\tilde{\Omega}}^{2}.$
(249)

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 $\breve{C}_{k}$, both independent of $e_{T}$, such
that

$\max_{x\in[0,l_{T}]}|\phi_{h}(x)-\phi(x)|\leq
C_{L,k}l_{T}^{k+2}|\phi|_{k+2,\infty,e_{T}}\leq\breve{C}_{k}l_{T}^{k+2}.$
(250)

Plugging (250) into (249) 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

$\begin{array}[]{l}\displaystyle\sum_{T\in{\mathcal{S}}_{0,h}}h_{T}^{-1}\|\tilde{u}\|_{0,\breve{e}_{T}}^{2}\leq\breve{C}_{L}\breve{C}_{k}^{2}h^{2k+3}\|\tilde{u}\|_{1,\infty,\tilde{\Omega}}^{2}\displaystyle\sum_{T\in{\mathcal{S}}_{0,h}}l_{T}\leq
C_{\tilde{\Omega}}^{2}length(\Gamma_{0})h^{2k+3}\|\tilde{u}\|_{3,\tilde{\Omega}}^{2}.\end{array}$
(251)

Finally, combining (251) and (246), we infer the existence of a constant
$\breve{C}_{\Gamma_{0}}$ depending only on $\Gamma_{0}$ and $\tilde{\Omega}$
such that the residual term $\breve{\mathcal{R}}_{h,k}(\tilde{u},{\bf q})$
fulfills

$|\breve{\mathcal{R}}_{h,k}(\tilde{u},{\bf
q})|\leq\breve{C}_{\Gamma_{0}}h^{k+3/2}\|\tilde{u}\|_{3,\tilde{\Omega}}[|{\bf
q}|]_{h}\;\forall{\bf q}\in{\bf Q}_{h,k}.$ (252)

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, ${\bf
p}\in[H^{k+1}(\Omega)]^{2}$, $u\in H^{\max[3,k+1]}(\Omega)$ and
$\nabla\cdot{\bf 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 the ingredients
exploited in Section 4, owing to the enhanced regularity of the underlying
auxiliary quantities.

##  7 Numerical experimentation

A comparative numerical study of the method addressed in this paper with the
classical (i.e. polygonal) approach was carried out in [26] for the cases
$k=0$ and $k=1$. Here we further do numerical experiments for $k=1$ taking a
non convex domain. Our aim is to highlight that optimal second order
convergence rates do apply to the case of a mixed Poisson problem with
Dirichlet conditions prescribed on a concave boundary portion $\Gamma_{0}$,
which confirms the predictions reported at the end of Subsection 7.1.
More specifically, we checked the performance of the formulation (10) for
$k=1$, referred to here as the P-G (for Petrov-Galerkin) $RT_{1}$ method, by
solving a test-problem in an annulus with inner radius $r_{i}=1/2$ and outer
radius $r_{e}=1$, for an exact solution given by
$u(x,y)=(r^{2}-2rr_{e}+2r_{i}r_{e}-r_{i}^{2})/2$, where
$r=\sqrt{x^{2}+y^{2}}$. This function satisfies $u=0$ on $\Gamma_{0}$, namely,
the circle given by $r=r_{i}$ and $\partial u/\partial r=0$ on $\Gamma_{1}$,
that is, along the circle given by $r=r_{e}$. We computed with a quasi-uniform
family of meshes for a quarter annulus, constructed for a quarter unit disk
with $2L^{2}$ triangles for $L=2^{m}$, by removing the $L^{2}/2$ triangles
fully contained in the disk with radius $1/2$. For simplicity we set $h=1/L$.
In the upper part of Table 1 we supply the approximation errors of $u$, ${\bf
p}$ and $\nabla\cdot{\bf p}$ measured in the norm of $L^{2}(\Omega_{h})$,
taking $m=2,3,4,5,6$. Although we did not study this kind of errors, in order
to further illustrate the properties of the Petrov-Galerkin approximation
(10), in the lower part of Table 1 we display the maximum errors of the
computed DOFs for both $u_{h}$ and ${\bf p}_{h}$, represented in the form of a
mesh-dependent semi-norm $|\cdot|_{0,\infty,h}$.

$h$ | $\longrightarrow$ | $1/4$ | $1/8$ | $1/16$ | $1/32$ | $1/64$
---|---|---|---|---|---|---
$\parallel u_{h}-u\parallel_{0,h}$ | $\longrightarrow$ | 0.28440E-2 | 0.70676E-3 | 0.17641E-3 | 0.44086E-4 | 0.11021E-4
$\parallel{\bf p}_{h}-{\bf p}\parallel_{0,h}$ | $\longrightarrow$ | 0.38543E-2 | 0.97914E-3 | 0.24598E-3 | 0.61577E-4 | 0.15400E-4
$\parallel\nabla\cdot({\bf p}_{h}-{\bf p})\parallel_{0,h}$ | $\longrightarrow$ | 0.82685E-2 | 0.21250E-2 | 0.53575E-3 | 0.13424E-3 | 0.33579E-4
|  |  |  |  |  |
$|u_{h}-u|_{0,\infty,h}$ | $\longrightarrow$ | 0.12974E-1 | 0.31851E-2 | 0.79104E-3 | 0.19721E-3 | 0.49252E-4
$|{\bf p}_{h}-{\bf p}|_{0,\infty,h}$ | $\longrightarrow$ | 0.66908E-2 | 0.17059E-2 | 0.44211E-3 | 0.11378E-3 | 0.29004E-4
Table 1: The P-G $RT_{1}$ method: errors for a test-problem in a non convex
curved domain

The predicted second order convergence in norm of our method in all the three
senses can be observed in Table 1. The same rate also applies to the computed
degrees of freedom of both $u_{h}$ and ${\bf p}_{h}$.

##  8 Final remarks and conclusions

We conclude this article with a few comments.

###  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}[]{l}\mbox{Given }f\in L^{2}(\Omega)\mbox{ find }({\bf
p};u)\in{\bf H}(div,\Omega)\times L^{2}(\Omega)\mbox{ such that}\\\
-(\nabla\cdot{\bf p},v)=(f,v)\;\forall v\in L^{2}(\Omega);\\\ ({\bf p},{\bf
q})+(u,\nabla\cdot{\bf q})=0\;\forall{\bf q}\in{\bf
H}(div,\Omega).\end{array}\right.$ (253)

Mimicking (10) and observing that we do not need ${\bf P}_{h,k}$ in this case,
(253) can be approximated by

$\left\\{\begin{array}[]{l}\mbox{Find }({\bf p}_{h};u_{h})\in{\bf
Q}_{h,k}\times V_{h,k}\mbox{ such that}\\\ -(\nabla\cdot{\bf
p}_{h},v)=(f,v)\;\forall v\in V_{h,k};\\\ ({\bf p}_{h},{\bf
q})+(u_{h},\nabla\cdot{\bf q})=0\;\forall{\bf q}\in{\bf
Q}_{h,k}.\end{array}\right.$ (254)

According to the predictions in Section 6, everything works fine for $k=0$ and
$k=1$, as long as $u\in H^{k+2}(\Omega)$. Otherwise stated, (254) yields
optimally converging approximations for $k\leq 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 ${\bf Q}_{h,k}\times V_{h,k}$ is
necessary in this case, i.e., the classical one given in [9]. However, if
$k>1$ the variational residual dominates the interpolation error, and hence
the optimal order in the sense of ${\bf H}(div,\Omega_{h}^{{}^{\prime}})\times
L^{2}(\Omega^{{}^{\prime}}_{h})$ can only be attained if we interpolate
$\Gamma_{T}$ by a polynomial arc of degree $k+1$ $\forall
T\in{\mathcal{S}}_{h}$ and integrate in the underlying curved element instead
of $T$, as explained in Subsection 6.3.

###  8.2 Summary and concluding remarks

This article is an outstanding application of the theory studied in [27] to
support the reliability analysis of numerical methods based on mixed
variational formulations with different product spaces for shape- and test-
fields. This tool is aimed at overcoming difficulties stemming from such a
symmetry breaking, as compared to the widespread theory of mixed formulations
with the same product spaces.

Akin to [24], [23] and [25], asymmetric Petrov-Galerkin formulations using
different spaces of shape- and test-functions are designed for the finite-
element solution of boundary value problems posed in smooth domains. The idea
behind this technique is the enforcement of essential conditions on the
boundary of the curved domain, while the test-functions remain attached to the
polytopal computational domain consisting of the union of the straight-edged
simplexes of a mesh fitting the true domain.

Although the application of our method to more complex problems is
straightforward, we studied it here in the particular case of the Poisson
equation in a curved plane domain with prescribed normal fluxes across the
boundary, in usual mixed formulation. Likewise previous work on this technique
listed above and in the bibliography, the proposed methodology brings about
enhanced accuracy, as compared to the formulation employing the same product
spaces of shape- and test-fields, defined upon an approximating polygonal
domain. For example, as reported in [26], in the case of the $RT_{1}$ method
it turns out that the mean-square errors of the flux variable are roughly
divided by two if our method is employed, while those of the multiplier remain
practically unchanged for both approaches.

An enriching topic for future work is the extension of our method to the
three-dimensional case in a context of a similar nature. This means
considering the Raviart-Thomas method for tetrahedrons.

## References

  * [1] R.A. Adams, Sobolev Spaces, Academic Press, N.Y., 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, Numerische Mathematik, 20 (1973), 179–192.
  * [4] F. Bertrand, S. Münzenmaier and G. Starke, First-order system least-squares on curved boundaries: higher-order Raviart–Thomas elements, SIAM J. Num. Analysis, 52-6 (2014), 3165–3180.
  * [5] F. Bertrand and G. Starke, Parametric Raviart-Thomas mixed methods on domains with curved surfaces, SIAM J. Numerical Analysis, 54-6 (2016), 3648–3667.
  * [6] E. K. Blum, Numerical Analysis and Computation: Theory and Practice, Addison Wesley, 1972.
  * [7] S.C. Brenner and L.R.Scott, The Mathematical Theory of Finite Element Methods, Texts Appl. Math. 15, Springer, 2008.
  * [8] F. Brezzi, On the existence, uniqueness and approximation of saddle-point problems arising from Lagrange multipliers, RAIRO Analyse Numérique, 8-2 (1974), 129-151.
  * [9] F. Brezzi and M. Fortin, Mixed and Hybrid Finite Element Methods, Springer-Verlag, 1991.
  * [10] 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., A.K. Aziz ed., pp. 409–474, Academic Press, 1972.
  * [11] P.G. Ciarlet, The Finite Element Method for Elliptic Problems. North Holland, 1978.
  * [12] P. G. Ciarlet, Three-dimensional elasticity, Vol. 1, Elsevier Science, 1994.
  * [13] J.A. Cuminato and V. Ruas, Unification of distance inequalities for linear variational problems, Computational and Applied Mathematics, 34 (2015), 1009-1033.
  * [14] E. Di Nezza, G. Palatucci, E. Valdinoci, Hitchhiker’s guide to the fractional Sobolev spaces, Bull. Sci. Math., 136 (2012), 521–573.
  * [15] G. Duvaut and J.-L. Lions,  Les Inéquations en Mécanique et en Physique, Dunod, 1972.
  * [16] A. Ern and J.L. Guermond,  Finite Elements I, Approximation and Interpolation, Springer, 2021.
  * [17] E.H. Georgoulis, Inverse-type estimates on $hp$-finite element spaces and applications, Mathematics of Computation 77-261 (2008), 201–219.
  * [18] V. Girault and P.A. Raviart, Finite Element Methods for Navier–Stokes Equations: Theory and Algorithms, Springer, 1986.
  * [19] P. Grisvard, Elliptic Problems in Non Smooth Domains, Pitman Publ. Inc, 1985.
  * [20] A. J. Lew and M. Negri, Optimal convergence of a discontinuous-Galerkin-based immersed boundary method, ESAIM M2AN, 45 (2011), 251–274.
  * [21] A. Quarteroni, R. Sacco and F. Saleri, Numerical Mathematics, Texts Appl. Math., Springer, 2007.
  * [22] P.-A. Raviart and J.-M. Thomas, Mixed Finite Element Methods for Second Order Elliptic Problems, Lecture Notes in Mathematics, Springer Verlag, 606: 292–315, 1977.
  * [23] V. Ruas, Optimal-rate nonconforming finite-element solution of Dirichlet problems in curved domains with straight-edged tetrahedra, IMA J. Num. Analysis, 41-2 (2021), 1368–1410.
  * [24] V. Ruas, Optimal Lagrange and Hermite finite elements for Dirichlet problems in curved domains with straight-edged triangles, ZAMM, 100-11 (2020), e-2019002962020.
  * [25] V. Ruas and M.A. Silva Ramos, Efficiency on non parametric finite elements for optimal-order enforcement of Dirichlet conditions on curvilinear boundaries, Journal of Computational and Applied Mathematics, 394 (October 2021; on-line), 113523.
  * [26] V. Ruas, Success and failure of attempts to improve the accuracy of Raviart–Thomas mixed finite elements in curved domains Examples and Counterexamples, Special Issue, Vol. 2 (November 2022; on-line), 100037.
  * [27] V. Ruas, Fine uniform stability constants and error bounds for asymmetric approximate saddle point problems, arXiv:2307.03742 [math.NA] (submitted paper), 2023.
  * [28] 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,  Numerische Mathematik, 45 (1984), 301–321.
  * [29] D. B. Stein, R. D. Guy, B. Thomases, Immersed boundary smooth extension: A high-order method for solving PDE on arbitrary smooth domains using Fourier spectral methods, Journal of Computational Physics, 304 (2016), 252–274.
  * [30] G. Strang and G. Fix,  An Analysis of the Finite Element Method, Prentice-Hall, 1973.
  * [31] J.L. Taylor, S. Kim and R.M. Brown, Comm. Partial Differential Eqns., 38-9 (2013), 1574-1600.
  * [32] R. Verfürth,  A Posteriori Error Estimation Techniques for Finite Element Methods, Oxford Science Publication, 2013.
  * [33] O.C. Zienkiewicz, The Finite Element Method in Engineering Science, McGraw-Hill, 1971.