\newsiamremark

remarkRemark \newsiamremarkhypothesisHypothesis \newsiamthmclaimClaim \headersAn a posteriori error estimate of the outer normal derivative using dual weightsS. Bertoluzza, E. Burman, and C. He

An a posteriori error estimate of the outer normal derivative using dual weights^†^†thanks: Submitted to the editors of SIAM Journal of Numerical Analysis. \fundingEB and CH were funded by the EPSRC grant EP/P01576X/1.

Silvia Bertoluzza Istituto di Matematica Applicata e Tecnologie Informatiche, CNR, Italy, (). [email protected] Erik Burman Department of Mathematics, University College London, UK, (). [email protected] Cuiyu He School of Mathematical and Statistical Sciences, University of Texas Rio Grand Valley, USA, (). [email protected]

Abstract

We derive a residual based a-posteriori error estimate for the outer normal flux of approximations to the diffusion problem with variable coefficient. By analyzing the solution of the adjoint problem, we show that error indicators in the bulk may be defined to be of higher order than those close to the boundary, which lead to more economic meshes. The theory is illustrated with some numerical examples.

keywords:

a posteriori error estimate; normal flux; dual weighted residual method

{AMS}

65M50, 65M60

Let ${\Omega}\subset\mathbb{R}^{d}$ , $d=2,3$ , be a polygonal/polyhedral domain, let $\Gamma=\partial{\Omega}$ denote its boundary and $\nu$ the outer unit normal. We consider the following diffusion problem

-\nabla\cdot{a}\nabla u=f,\mbox{ in }\Omega,

with non homogeneous Dirichlet boundary conditions, $u=g$ on $\Gamma$ . The outer normal flux $\nu\cdot({a}\nabla u)$ is an important quantity in many applications. It is of importance for instance when a heat flux or an electric field on the boundary of the domain needs to be approximated, or in fluid mechanics for the fluid forces [1, 29, 36, 20]. For boundary control problems, an accurate approximation of the normal flux on the boundary also plays a critical role [2, 3]. Recently there has been a number of works estimating the error for the outer normal flux in the a priori sense. We refer to [32, 35].

From the computational perspective it is appealing to apply adaptive methods that concentrate degrees of freedom where they are most needed to achieve a certain accuracy. In particular, for the normal flux on the boundary, we expect perturbations in the bulk of the domain to be less significant than those close to the boundary. This is proved in [14] where local a priori error estimates were given for the error in the outer normal flux. In particular, the error on the flux quantity was shown to depend on the $H^{1}$ -error in a tubular neighborhood of the boundary and a global term that measures the global error in a weak norm. Similar results using boundary concentrated meshes were obtained more recently in [38], where the application to a Dirichlet boundary control problem was studied. A consequence of the localization property underlying the above a priori error estimates is that a standard energy norm estimate is unlikely to have optimal performance when approximating the normal flux, since it does not account for the relative independence of the goal quantity on perturbations in the bulk. It is however not straightforward to ensure accuracy of the boundary flux using a priori refinement in the boundary region alone, since geometric singularities or rough data nevertheless have to be taken into account.

The objective of the present work is to derive a residual based a posteriori error estimate for the outer normal flux that exploits the localization property. In particular, we add some mesh dependent weight in front of the classical residual based error estimator, and the weights greatly depend on the distance to the boundary. More precisely, the domain is implicitly divided into two zones, a tubular neighborhood around the boundary and an interior, bulk zone. For elements in the latter, the residual estimator is multiplied with the mesh diameter to a higher power than in the boundary region, hence giving it relative smaller weight. To get a precise quantification of the size of the weight we consider an adjoint problem. Thanks to suitable weighted estimates we determine the rate of the decrease of the adjoint solution and its derivatives with increasing distance to the boundary. This then helps provide bounds on the dual weights in the a posteriori error estimate that allow us to decompose the domain in a bulk and a boundary subdomain with associated error indicators.

The use of adjoint equations for the derivation of a posteriori error estimates in weak norms was first proposed by Eriksson and Johnson in [21], in the case of $L^{2}$ -norm bounds. These ideas were generalized to the approximation of fluxes and fluid forces using the Dual Weighted Residual a posteriori error estimation approach (see for instance [9, 10, 27, 8, 11, 39]). In these approaches, the dual solution was approximated, typically focussing on linear functionals of the error. There has recently been an increased interest in the convergence and optimality of goal oriented adaptive methods [5, 31, 24, 30, 7, 33, 6]. With this work we show that when the target quantity of the computation is the outward normal flux, a detailed analysis of the adjoint equation can lead to a posteriori bounds that perform better than the standard energy estimate, but without the need of solving the dual problem, numerically. Recall in this context that, when the target quantity driving the adaptive procedure is a norm of the error, the computation of the solution of the adjoint problem is complicated by the fact that the right hand side depends on the error itself, and is therefore not directly available, contrary to what happens for instance when the target quantity is a given known functional of the solution, such as the value at a point or the integral over a line.

Herein we only consider the standard finite element setting where the domain is meshed with a conforming triangulation. However, the arguments generalize in a straightforward manner to a posteriori error estimates for fictitious domain methods where elements are cut [17]. To extend the method to adaptive standard fictitious domain methods [28] in the spirit of [13], or domain decomposition methods, some more subtle arguments are needed. Indeed in such situations, the boundary divides the computational domain in two (or more) subdomains, thus requiring an analysis of the adjoint solution, similar to the one in this paper, for each subdomain and accounting for all boundaries and interfaces of the problem. This is the topic of a forthcoming paper.

An outline of the paper is as follows. First we introduce the weak formulation of our model problem and the associated finite element method in section 1. In section 2 we derive the a posteriori error estimate. Then we show in section 3 how to apply the results to some known stabilized methods, such as the Barbosa-Hughes methods and Nitsche’s method. Finally, we illustrate the theory with some numerical examples in section 4.

1 The Lagrange multiplier formulation of the Dirichlet Problem

For $g\in H^{1/2}(\Gamma)$ and $f\in L^{2}(\Omega)$ given, we consider the problem of finding $u\in H^{1}({\Omega})$ , ${\lambda}\in H^{-1/2}(\Gamma)$ such that for all $v\in H^{1}({\Omega})$ , $\mu\in H^{-1/2}(\Gamma)$

(1)

\displaystyle\int_{\Omega}a\,\nabla u\cdot\nabla v-\int_{\Gamma}{\lambda}v=\int_{\Omega}fv,\qquad\qquad\int_{\Gamma}u\mu=\int_{\Gamma}g\mu.

where $a\in C^{\infty}(\bar{\Omega})$ is the diffusion coefficient, which for the sake of simplicity we assume to be scalar, satisfying $0<\alpha\leq a\leq M$ for some constants $\alpha$ and $M$ . We consider a Galerkin discretization of such problem. More precisely, letting $V_{h}\subset H^{1}({\Omega})$ , ${\Lambda}_{h}\subset H^{-1/2}(\Gamma)$ be finite element spaces defined on a shape regular triangulation ${\mathcal{T}}_{h}$ . We look for $u_{h}\in V_{h}$ , ${\lambda}_{h}\in{\Lambda}_{h}$ such that for all $v_{h}\in V_{h}$ , $\mu_{h}\in{\Lambda}_{h}$

(2)

\int_{\Omega}a\nabla u_{h}\cdot\nabla v_{h}-\int_{\Gamma}{\lambda}_{h}v_{h}=\int_{\Omega}fv_{h},\qquad\qquad\int_{\Gamma}u_{h}\mu_{h}=\int_{\Gamma}g\mu_{h}.

We assume that $V_{h}$ contains the space of continuous piecewise polynomials of order $k$ $(k\geq 0)$ on ${\mathcal{T}}_{h}$ , which we denote by $\widecheck{V}_{h}$ , and that $\Lambda_{h}$ contains a subspace $\widecheck{\Lambda}_{h}$ which is either the space of piecewise constants, or the space of continuous piecewise linears on the mesh induced on $\Gamma$ by ${\mathcal{T}}_{h}$ .

Restricting the test functions in Eq. 1 to the discrete spaces and taking the difference of Eq. 1 and Eq. 2 we see that the following Galerkin orthogonality holds: for all $v_{h}\in V_{h}$ , $\mu_{h}\in{\Lambda}_{h}$

(3)

\int_{\Omega}a\,\nabla(u-u_{h})\cdot\nabla v_{h}-\int_{\Gamma}(\lambda-{\lambda}_{h})v_{h}=0,\qquad\qquad\int_{\Gamma}(u-u_{h})\mu_{h}=0.

Observe that in the above we are as general as possible in the definition of the two spaces. We do not even need to assume that the spaces satisfy the inf-sup condition required for the stability of Eq. 2. This of course does not mean that the method is stable without it, only that the a posteriori error estimate will measure the computational error independently of the stability properties of the pair $V_{h}\times\Lambda_{h}$ . An example of spaces that may be used in the framework are

(4)

V_{h}=\{u\in H^{1}(\Omega):\ u|_{T}\in\mathbb{P}_{k}(T),\ \forall T\in{\mathcal{T}}_{h}\},

and, for $k^{\prime}\geq 0$

(5)

\Lambda_{h}=\{\lambda\in L^{2}(\Gamma):\ u|_{F}\in\mathbb{P}_{k^{\prime}}(F),\ \forall F\in{\mathcal{T}}_{h}|_{\Gamma}\},

or, for $k^{\prime}\geq 1$

(6)

\Lambda_{h}=\{\lambda\in C^{0}(\Gamma):\ u|_{F}\in\mathbb{P}_{k^{\prime}}(F),\ \forall F\in{\mathcal{T}}_{h}|_{\Gamma}\}.

Also variants of the spaces Eq. 5 and Eq. 6 with local conforming enrichment on the boundary to satisfy the inf-sup condition are valid [16].

Remark 1.1.

We point out that, for $k^{\prime}=k$ , the choice Eq. 6 for the multiplier space, coupled with the choice Eq. 4 for the approximation of the primal unknown (i.e., choosing $\Lambda_{h}=V_{h}|_{\Gamma}$ ) yields a stable discretization of Problem Eq. 1, equivalent to strongly imposing the Dirichlet boundary condition $u_{h}=\pi_{h}g$ , where $\pi_{h}:L^{2}(\Gamma)\to V_{h}|_{\Gamma}$ is the $L^{2}(\Gamma)$ orthogonal projection. Then, using $\lambda_{h}$ as an approximation to the normal flux is equivalent to compute the latter by post-processing with a variational approach as proposed, for instance, in [38]. Remark that, when the domain has corners, this method will not have, in general, optimal approximation for the multiplier, and it should be modified following the strategy used in the mortar method (see [12]), where discontinuity is allowed at the corners, with $k^{\prime}=k-1$ for those elements on the boundary mesh ${\mathcal{T}}_{h}|_{\Gamma}$ which are adjacent to the corners, and $k^{\prime}=k$ for the remaining elements. Observe, however, that also for the suboptimal choice Eq. 6, the estimator we are going to present, remains valid.

2 A posteriori error estimates

The a posteriori error estimate is derived in three steps. We first derive an error representation using the adjoint problem. We then derive the local bounds for the adjoint solution and, finally, we obtain the weighted residual estimates. In what follows we will use the notation $A\lesssim B$ to indicate that $A\leq cB$ for some positive constant $c$ independent of mesh size parameters such as element diameters and/or face diameters or edge lengths. $A\simeq B$ will stand for $A\lesssim B\lesssim A$ .

2.1 Error representation using duality

We let

A:(H^{1}(\Omega)\times H^{-1/2}(\Gamma))\times(H^{1}(\Omega)\times H^{-1/2}(\Gamma))\to\mathbb{R}

be defined by

(7)

A(w,\eta;v,\zeta)=\int_{\Omega}{a\,}\nabla w\cdot\nabla v-\int_{\Gamma}\eta v+\int_{\Gamma}w\zeta.

Let $(u,\lambda)\in H^{1}(\Omega)\times H^{-1/2}(\Gamma)$ be the solution of Eq. 1 and let $(u_{h},\lambda_{h})\in V_{h}\times\Lambda_{h}$ satisfy Eq. 2. Set $e=u-u_{h}$ and $\delta=\lambda-\lambda_{h}$ . We define $L:H^{-1/2}(\Gamma)\to\mathbb{R}$ as

L(\xi):=\|\delta\|^{-1}_{-1/2,\Gamma}(\delta,\xi)_{-1/2,\Gamma},\qquad\text{so that }\qquad L(\delta)=\|\delta\|_{-1/2,\Gamma},

where $(\cdot,\cdot)_{-1/2,\Gamma}$ is the scalar product for the space $H^{-1/2}(\Gamma)$ , whose precise expression is provided later in Eq. 11, and where $\|\cdot\|_{-1/2,\Gamma}$ is the corresponding norm. Define $(z,\zeta)\in{\cal V}=H^{1}(\Omega)\times H^{-1/2}(\Gamma)$ as the solution of

(8)

A(w,\eta;z,\zeta)=L(\eta),\qquad\forall\ (w,\eta)\in{\cal V}.

Remark that the right hand side functional $L$ depends on the unknown error $\delta$ , so that it is not possible to compute $z,\zeta$ , even only approximately. It is, however, easy to see that $|L(\xi)|\leq\|\xi\|_{-1/2,\Gamma}$ , and then the operator $L$ has unitary norm. Therefore, by the stability of Eq. 8, we have

(9)

\|z\|_{1,\Omega}\lesssim 1,\qquad\|\zeta\|_{-1/2,\Gamma}\lesssim 1.

Let ${{\mathcal{F}_{h}^{i}}}$ and ${{\mathcal{F}_{h}^{b}}}$ respectively denote the set of interior and boundary $(d-1)$ -dimensional facets of the triangulation ${\mathcal{T}}_{h}$ and, for an element $T\in{\mathcal{T}}_{h}$ , let $\nu_{T}$ denote the outer unit normal to $\partial T$ . On a $(d-1)$ -dimensional facet $F=\partial T^{+}\cap\partial T^{-}$ we define the jump of the normal flux by $\Lbrack{a\,}\partial_{\nu}u_{h}\Rbrack={a\,}\nabla u_{h}^{+}\cdot\nu_{T^{+}}+{a\,}\nabla u_{h}^{-}\cdot\nu_{T^{-}}$ .

Proposition 2.1.

(Error representation) Let $\delta=\lambda-\lambda_{h}$ and let $z,\zeta$ be the solution of Eq. 8. Then it holds that for any $z_{h}\in V_{h}$ and $\zeta_{h}\in\Lambda_{h}$

(10)

\begin{split}\|\delta\|_{-1/2,\Gamma}=&\sum_{T\in{\mathcal{T}}_{h}}\int_{T}(f+{\nabla\cdot a\,\nabla}u_{h})(z-z_{h})-\sum_{F\in{{\mathcal{F}_{h}^{i}}}}\int_{F}\Lbrack{a\,}\partial_{\nu}u_{h}\Rbrack(z-z_{h})\\ &+\sum_{F\in{{\mathcal{F}_{h}^{b}}}}\int_{F}(\lambda_{h}-{a\,}\partial_{\nu}u_{h})(z-z_{h})+\int_{\Gamma}(g-u_{h})(\zeta-\zeta_{h}).\end{split}

Proof 2.2.

Taking $w=e$ and $\eta=\delta$ in Eq. 8 we have

\|\delta\|_{-1/2,\Gamma}=L(\delta)=A(e,\delta;z,\zeta).

Now, for $z_{h}\in V_{h}$ , $\zeta_{h}\in\Lambda_{h}$ arbitrary, thanks to Galerkin orthogonality Eq. 3 we can write:

\|\delta\|_{-1/2,\Gamma}=A(e,\delta;z-z_{h},\zeta-\zeta_{h})=I+II

with

I=\int_{\Omega}a\nabla e\cdot\nabla(z-z_{h}),\;II=-\int_{\Gamma}(a\partial_{\nu}u-\lambda_{h})(z-z_{h})+\int_{\Gamma}(g-u_{h})(\zeta-\zeta_{h}).

For the term $I$ we obtain using Green’s theorem

\begin{split}I=&\int_{\Omega}a\nabla e\cdot\nabla(z-z_{h})=\sum_{T\in{\mathcal{T}}_{h}}\int_{T}aa\nabla e\cdot\nabla(z-z_{h})\\ =&\sum_{T\in{\mathcal{T}}_{h}}\left(\int_{T}(f+{\nabla\cdot a\,\nabla}u_{h})(z-z_{h})+\int_{\partial T}a\nabla(u-u_{h})\cdot\nu_{T}(z-z_{h})\right)\\ =&\sum_{T\in{\mathcal{T}}_{h}}\int_{T}(f+{\nabla\cdot a\,\nabla}u_{h})(z-z_{h})-\sum_{F\in{{\mathcal{F}_{h}^{i}}}}\int_{F}\Lbrack a\partial_{\nu}u_{h}\Rbrack(z-z_{h})\\ &+\sum_{F\in{{\mathcal{F}_{h}^{b}}}}\int_{F}(a\partial_{\nu}u-a\partial_{\nu}u_{h})(z-z_{h}).\end{split}

Combining all yields Eq. 10. This completes the proof of the proposition.

2.1.1 Some observations on the operator $L$

We start by observing that taking $v_{h}=1$ in Eq. 3 implies $\int_{\Gamma}\delta=0$ . Then we have

\|\delta\|_{-1/2,\Gamma}=\sup_{\phi\in H^{1/2}(\Gamma)}\frac{\int_{\Gamma}\delta\phi}{\|\phi\|_{1/2,\Gamma}}\simeq\sup_{{\phi\in H^{1/2}(\Gamma)}\atop{\int_{\Gamma}\phi=0}}\frac{\int_{\Gamma}\delta\phi}{|\phi|_{1/2,\Gamma}}.

On the space $H^{1/2}_{\circ}(\Gamma)=\{\phi\in H^{1/2}:\ \int_{\Gamma}\phi=0\}$ of zero average functions in $H^{1/2}(\Gamma)$ , we can define a scalar product and a norm, equivalent to the standard $H^{1/2}$ scalar product and norm, as

(\phi,\psi)_{1/2,\Gamma}=\int_{\Omega}\nabla\phi^{\mathcal{H}}\cdot\nabla\psi^{\mathcal{H}},\qquad|\phi|_{1/2,\Gamma}:=|\phi^{\mathcal{H}}|_{1,\Omega},

where $\phi^{\mathcal{H}}\in H^{1}(\Omega)$ denotes the harmonic lifting of $\phi$ . We then let $\|\cdot\|_{-1/2,\Gamma}$ be defined by duality with respect to the above norm. We now let $\mathfrak{R}:(H^{1/2}_{\circ}(\Gamma))^{\prime}\to H^{1/2}_{\circ}(\Gamma)$ denote the Riesz isomorphism, which, we recall, is defined as the solution of

(\mathfrak{R}\lambda,\phi)_{1/2,\Gamma}=\int_{\Gamma}\lambda\phi\quad\forall\,\phi\in H^{1/2}_{\circ}(\Gamma).

We recall that, as $\mathfrak{R}$ is an isomorphism, we also have that

(11)

{(\lambda,\mu)_{-1/2,\Gamma}=(\mathfrak{R}\lambda,\mathfrak{R}\mu)_{1/2,\Gamma}}=\int_{\Omega}\nabla(\mathfrak{R}\lambda)^{\mathcal{H}}\cdot\nabla(\mathfrak{R}\mu)^{\mathcal{H}}.

It is now easy to check that, if $\mu\in L^{2}(\Gamma)$ satisfies $\int_{\Gamma}\mu=0$ , then $(\mathfrak{R}\mu)^{\mathcal{H}}$ is the unique solution of

(12)

-\Delta(\mathfrak{R}\mu)^{\mathcal{H}}=0\text{ in }\Omega,\qquad\int_{\Gamma}(\mathfrak{R}\mu)^{\mathcal{H}}=0,\qquad\partial(\mathfrak{R}\mu)^{\mathcal{H}}/\partial\nu=\mu.

Indeed for any function $v\in H^{1}(\Omega)$ , there is a unique decomposition $v=\bar{v}+v_{1}+v_{0}$ such that $\bar{v}=|\Gamma|^{-1}\int_{\Gamma}v$ , $v_{0}\in H^{1/2}_{\circ}(\Gamma)$ is the harmonic extension of $v-\bar{v}$ , and $v_{1}\in H_{0}^{1}(\Omega)$ satisfies $\triangle v_{1}=\triangle v$ . Then we have that for any $v\in H^{1}(\Omega)$

(13)

\begin{split}&\int_{\Omega}\nabla(\mathfrak{R}\mu)^{\mathcal{H}}\cdot\nabla v=\int_{\Omega}\nabla(\mathfrak{R}\mu)^{\mathcal{H}}\cdot\nabla v_{0}=(\mathfrak{R}\mu,v_{0})_{1/2}=\int_{\Gamma}\mu v_{0}=\int_{\Gamma}\mu v,\end{split}

which is the weak form of equation Eq. 12.

2.2 Local estimates for the adjoint solution $z$

We observe that $z$ is the solution of the following problem.

\int_{\Omega}a\nabla w\cdot\nabla z+\int_{\Gamma}w\zeta=0,\qquad-\int_{\Gamma}\eta z=\|\delta\|_{-1/2,\Gamma}^{-1}(\delta,\eta)_{-1/2,\Gamma}=|\mathfrak{R}\delta|^{-1}_{1/2,\Gamma}\int_{\Gamma}\eta\,\mathfrak{R}\delta.

This rewrites as

-{\nabla\cdot a\,\nabla}z=0\text{ in }\Omega,\qquad z=-|\mathfrak{R}\delta|^{-1}_{1/2,\Gamma}{\mathfrak{R}\delta}{}\text{ on }\Gamma.

The following Lemma, whose proof we include for the sake of completeness, was proven in [34].

Lemma 2.3.

Let $d_{\Gamma}(x)$ denote the distance of $x$ from $\Gamma$ and let $w\in H^{1}(\Omega)$ satisfy ${\nabla\cdot a\,\nabla}w=0$ in $\Omega$ . Then, for all $p\geq 0$ it holds that

(14)

\|d_{\Gamma}^{p+1}\nabla^{p+2}w\|_{0,\Omega}\lesssim|w|_{1,\Omega}.

Proof 2.4.

We start by proving a local bound. Let $B_{R}$ and $B_{cR}$ , $0<c<1$ be two concentric balls of radius respectively $R$ and $cR$ , and assume that $w\in H^{1}(B_{R})$ satisfies ${\nabla\cdot a\,\nabla}w=0$ in $B_{R}$ . Then, we claim that for all $p\geq 0$ it holds that

(15)

\|\nabla^{p+2}w\|_{0,B_{cR}}\lesssim R^{-p-1}\|\nabla w\|_{0,B_{R}}+R^{-p-2}\|w\|_{0,B_{R}},

where the implicit constant in the inequality depends on $c$ . We start by provingEq. 15 for $R=1$ . We prove it by induction on $p$ . For $p=0$ , this is a consequence of [26, Theorem 8.8]. Let us now assume that the result is true for all $p\leq n-1$ and prove it for $p=n$ . We let $c^{\prime}=1-(1-c)/2=c/2+1/2$ , and let $\omega_{c}\in C^{\infty}_{0}(B_{c^{\prime}})$ , $\omega_{c}\geq 0$ , $\omega_{c}=1$ in $B_{c}$ . We have

{\nabla\cdot a\,\nabla}(\omega_{c}w)=2a\nabla w\cdot\nabla\omega_{c}+aw\Delta\omega_{c}+w\nabla a\cdot\nabla\omega_{c},\quad\text{ in }B_{c^{\prime}}\qquad\omega_{c}w=0,\quad\text{ on }\partial B_{c^{\prime}}.

Using standard results on the smoothness of the solution of elliptic equations (see [26]), by the induction assumption we have that

\begin{split}\|\nabla^{n+2}w\|_{0,B_{c}}&\leq\|\nabla^{n+2}(\omega_{c}w)\|_{0,B_{c^{\prime}}}\lesssim\|2a\nabla w\cdot\nabla\omega_{c}+aw\Delta\omega_{c}+w\nabla a\cdot\nabla\omega_{c}\|_{n,B_{c^{\prime}}}\\ &\lesssim\|\nabla w\|_{n,B_{c^{\prime}}}+\|w\|_{n,B_{c^{\prime}}}\lesssim\|w\|_{1,B_{1}},\end{split}

which proves our claim for $R=1$ . By rescaling we immediately obtain Eq. 15.

Let us now prove Eq. 14. We consider a covering of $\Omega$ , consisting of a countable collection of balls $B_{i}=B_{r_{i}}(x_{i})\subset\Omega$ , of center $x_{i}$ and radius $r_{i}$ , with $r_{i}=\tilde{c}d_{\Gamma}(x_{i})$ for some fixed $0<\tilde{c}<1$ , such that

1.

there exist $N\in\mathbb{N}$ such that all $x\in\Omega$ belong to at most $N$ balls $B_{i}$ ;
2.

for some $0<c<1$ independent of $i$ , letting ${\widetilde{B}_{i}}\subset\subset B_{i}$ denote the ball of center $x_{i}$ and radius $cr_{i}$ , it holds that $\Omega\subseteq\cup_{i}{\widetilde{B}_{i}}$ ,

(by the Besicovitch covering Theorem, such a collection exists). We observe that the relation between the radius of the balls in our covering and the distance of the centers from the boundary of the domain implies that for all $i$ , $x\in{\widetilde{B}_{i}}$ implies $d_{\Gamma}(x)\simeq r_{i}$ . Then, letting $w_{i}=|B_{i}|^{-1}\int_{B_{i}}w$ denote the average of $w$ in $B_{i}$ , using Eq. 15 and a Poincaré inequality, we can write

\begin{split}\|d_{\Gamma}^{p+1}\nabla^{p+2}w\|_{0,\Omega}^{2}&\leq\sum_{i}\|d_{\Gamma}^{p+1}\nabla^{p+2}w\|_{0,{\widetilde{B}_{i}}}^{2}\lesssim\sum_{i}r_{i}^{2(p+1)}\|\nabla^{p+2}w\|_{0,{\widetilde{B}_{i}}}^{2}\\ &\lesssim\sum_{i}r_{i}^{2(p+1)}\|\nabla^{p+2}(w-w_{i})\|_{0,{\widetilde{B}_{i}}}^{2}\lesssim\sum_{i}(|w-w_{i}|^{2}_{1,B_{i}}+r_{i}^{-2}\|w-w_{i}\|_{0,B_{i}}^{2})\\ &\lesssim\sum_{i}|w|_{1,B_{i}}^{2}\lesssim|w|_{1,\Omega}^{2},\end{split}

which concludes the proof.

2.3 The a posteriori error estimator

Using the error representation of Proposition 2.1 and the local bounds for the adjoint solution stated in Eq. 14, we will now derive the a posteriori error estimation. Comparing to the classical residual based error indicator, our local error indicators for each element/facet are additionally multiplied by local dual weights depending on the distance from the element/facet to the boundary. Let us first introduce some notations that will be useful for the bounds.

We let $h_{T}$ (resp. $h_{F}$ ) denote the diameter of an element $T$ (resp. of a $(d-1)$ -dimensional facet $F$ ) in ${\mathcal{T}}_{h}$ . For a given element $T\in{\mathcal{T}}_{h}$ , ${\Delta_{T}}$ denotes the patch of elements that have at least a vertex in common with $T$ . The distance of an element $T$ to the boundary will be measured using $\rho_{T}=\underset{x\in\Delta_{T}}{\min}d_{\Gamma}(x)$ . That is the shortest distance from the associated patch to the boundary.

We now let $\widehat{\Pi}_{h}:H^{1}(\Omega)\to\widecheck{V}_{h}$ denote the Scott-Zhang projector, introduced in [40]. We recall that, for $1\leq m\leq k+1$ it holds that

(16)

\displaystyle\|z-\widehat{\Pi}_{h}z\|_{0,T}+h_{T}|z-\widehat{\Pi}_{h}z|_{1,T}\lesssim h^{m}_{T}|z|_{m,{\Delta_{T}}}.

Using this bound for $m=1$ and $m=k+1$ we have the following local interpolation bounds for the adjoint solution.

Lemma 2.5.

Let $z_{h}=\widehat{\Pi}_{h}z$ , then we have the following two bounds

	$\displaystyle\\|z-z_{h}\\|_{0,T}+h_{T}\|z-z_{h}\|_{1,T}\leq C_{1}h_{T}\|z\|_{1,{\Delta_{T}}},$
	$\displaystyle{\\|z-z_{h}\\|_{0,T}+h_{T}\|z-z_{h}\|_{1,T}\leq C_{2}h_{T}^{k+1}\rho_{T}^{-k}\\|d_{\Gamma}^{k}\,\nabla^{k+1}z\\|_{0,{\Delta_{T}}}}$

The constants $C_{1}$ and $C_{2}$ depend on the shape regularity of the mesh

Proof 2.6.

The first inequality is immediate by Eq. 16 with $m=1$ . The second inequality trivially holds for $T$ with $\Delta_{T}$ adjacent to the boundary (for which $\rho_{T}^{-k}=\infty$ ). For the elements for which $\Delta_{T}$ is interior to $\Omega$ , it follows by first applying Eq. 16 with $m=k+1$ , then multiplying and dividing by $d_{\Gamma}^{k}$ , and finally bounding $d_{\Gamma}^{-k}\leq\rho_{T}^{-k}$ :

\|z-z_{h}\|_{0,T}+h_{T}|z-z_{h}|_{1,T}\lesssim h^{k+1}_{T}|z|_{k+1,{\Delta_{T}}}\lesssim h^{k+1}_{T}\rho_{T}^{-k}|d_{\Gamma}^{k}z|_{k+1,{\Delta_{T}}}.

Let us at first assume that we have $g\in H^{1}(\Gamma)$ . Under such an assumption we have the following theorem.

Theorem 2.7.

Define the following local residuals:

(17)

\begin{split}&{\mathbb{r}}(T)=h_{T}\|f+{\nabla\cdot a\,\nabla}u_{h}\|_{0,T},\quad\forall T\in\mathcal{T}_{h},\\ &\mathbf{r}_{0}(F)=h_{F}^{1/2}\|\Lbrack a\partial_{\nu}u_{h}\Rbrack\|_{0,F},\quad\forall F\in{{\mathcal{F}_{h}^{i}}},\\ &\mathbf{r}_{1}(F)=h_{F}^{1/2}\|\lambda_{h}-a\partial_{\nu}u_{h}\|_{0,F},\quad\forall F\in{{\mathcal{F}_{h}^{b}}},\\ &\mathbf{r}_{2}(F)=h_{F}^{1/2}|g-u_{h}|_{1,F},\quad\forall F\in{{\mathcal{F}_{h}^{b}}}.\end{split}

Then we have

(18)

\|\lambda-\lambda_{h}\|_{-1/2,\Gamma}\lesssim\sqrt{\sum_{T\in{\mathcal{T}}_{h}}\varsigma_{T}^{2}|{\mathbb{r}}(T)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{i}}}}\varsigma_{F}^{2}|\mathbf{r}_{0}(F)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{b}}}}\left(|\mathbf{r}_{1}(F)|^{2}+|\mathbf{r}_{2}(F)|^{2}\right)},

where the element and facet weights $\varsigma_{T}$ and $\varsigma_{F}$ are defined by

(19)

\varsigma_{T}=\min\{C_{1},C_{2}h^{k}_{T}{\rho_{T}^{{-k}}}\},\quad{\varsigma_{F}=\min\{\varsigma_{T},\varsigma_{T^{\prime}}\}},\mbox{ with }F=T\cap T^{\prime}.

Proof 2.8.

Let us start by splitting ${\mathcal{T}}_{h}$ as the union of two disjoint sets

{\mathcal{T}}_{h}^{1}=\Big{\{}T\in{\mathcal{T}}_{h}:{C_{1}|z|_{1,{\Delta_{T}}}\leq C_{2}h_{T}^{k}\rho_{T}^{-k}\|d_{\Gamma}^{k}\nabla^{k+1}z\|_{0,{\Delta_{T}}}}\Big{\}},\qquad{\mathcal{T}}_{h}^{2}={\mathcal{T}}_{h}\setminus{\mathcal{T}}_{h}^{1}.

Setting $z_{h}=\widehat{\Pi}_{h}z$ and $\zeta_{h}=0$ in the error representation of Proposition 2.1, we have

\begin{split}\|\delta\|_{-1/2,\Gamma}=&\sum_{T\in{\mathcal{T}}_{h}}\int_{T}(f+{\nabla\cdot a\,\nabla}u_{h})(z-\widehat{\Pi}_{h}z)-\sum_{F\in{{\mathcal{F}_{h}^{i}}}}\int_{F}\Lbrack a\partial_{\nu}u_{h}\Rbrack(z-\widehat{\Pi}_{h}z)\\ &+\sum_{F\in{{\mathcal{F}_{h}^{b}}}}\int_{F}(\lambda_{h}-a\partial_{\nu}u_{h})(z-\widehat{\Pi}_{h}z)+\int_{\Gamma}(g-u_{h})\zeta.\end{split}

Observe that Lemma 2.5 gives us two error estimates for $\|z-\widehat{\Pi}_{h}z\|_{0,T}$ , and, depending on whether $T\in{\mathcal{T}}_{h}^{1}$ or $T\in{\mathcal{T}}_{h}^{2}$ , we apply the best possible estimate. This yields

\begin{split}&\sum_{T\in{\mathcal{T}}_{h}}\int_{T}(f+{\nabla\cdot a\,\nabla}u_{h})(z-\widehat{\Pi}_{h}z)\\ \lesssim&\sum_{T\in{{\mathcal{T}}_{h}^{1}}}\|f+{\nabla\cdot a\,\nabla}u_{h}\|_{0,T}C_{1}h_{T}|z|_{1,{\Delta_{T}}}+\sum_{T\in{{\mathcal{T}}_{h}^{2}}}\|f+{\nabla\cdot a\,\nabla}u_{h}\|_{0,T}C_{2}{h_{T}^{k}\rho_{T}^{-k}\|d_{\Gamma}^{k}\nabla^{k+1}z\|_{0,{\Delta_{T}}}}\\ =&\sum_{T\in{\mathcal{T}}_{h}^{1}}\varsigma_{T}{\mathbb{r}}(T)|z|_{1,{\Delta_{T}}}+\sum_{T\in{\mathcal{T}}_{h}^{2}}\varsigma_{T}{\mathbb{r}}(T){\|d_{\Gamma}^{k}\nabla^{k+1}z\|_{0,{\Delta_{T}}}}\\ \leq&\sqrt{\sum_{T\in{\mathcal{T}}_{h}}\varsigma_{T}^{2}|{\mathbb{r}}(T)|^{2}}\sqrt{\sum_{T\in{\mathcal{T}}_{h}^{1}}|z|_{1,{\Delta_{T}}}^{2}+{\sum_{T\in{\mathcal{T}}_{h}^{2}}\|d_{\Gamma}^{k}\nabla^{k+1}z\|_{0,{\Delta_{T}}}^{2}}}.\end{split}

Applying Lemma 2.3 we have

\sum_{T\in{\mathcal{T}}_{h}^{1}}|z|_{1,{\Delta_{T}}}^{2}+{\sum_{T\in{\mathcal{T}}_{h}^{2}}\|d_{\Gamma}^{k}\nabla^{k+1}z\|_{0,{\Delta_{T}}}^{2}}\lesssim\|z\|_{1,\Omega}^{2}+{\|d_{\Gamma}^{k}\nabla^{k+1}z\|_{0,\Omega}^{2}\lesssim\|z\|_{1,\Omega}}\lesssim 1,

so that

(20)

\sum_{T\in{\mathcal{T}}_{h}}\int_{T}(f+{\nabla\cdot a\,\nabla}u_{h})(z-\widehat{\Pi}_{h}z)\lesssim\sqrt{\sum_{T\in{\mathcal{T}}_{h}}\varsigma_{T}^{2}|{\mathbb{r}}(T)|^{2}},

where the constant in the inequality depends on $\Omega$ and $a$ .

A similar argument can be applied for interior facets. Letting $F\in{{\mathcal{F}_{h}^{i}}}$ , $F\subset\partial T$ , the standard bound holds

(21)

\begin{split}&\int_{F}\Lbrack a\partial_{\nu}u_{h}\Rbrack(z-\widehat{\Pi}_{h}z)\leq\|\Lbrack a\partial_{\nu}u_{h}\Rbrack\|_{0,F}\|z-\widehat{\Pi}_{h}z\|_{0,F}\\ \lesssim&\|\Lbrack a\partial_{\nu}u_{h}\Rbrack\|_{0,F}\left(h_{T}^{-1/2}\|z-\widehat{\Pi}_{h}z\|_{0,T}+h_{T}^{1/2}|z-\widehat{\Pi}_{h}z|_{1,T}\right)\leq\|\Lbrack a\partial_{\nu}u_{h}\Rbrack\|_{0,F}h_{T}^{1/2}C_{1}|z|_{1,{\Delta_{T}}},\end{split}

as well as the enhanced bound

(22)

\int_{F}\Lbrack a\partial_{\nu}u_{h}\Rbrack(z-\widehat{\Pi}_{h}z)\leq\|\Lbrack a\partial_{\nu}u_{h}\Rbrack\|_{0,F}C_{2}h_{T}^{k+1/2}{\rho_{T}^{-k}\|d_{\Gamma}^{k}\nabla^{k+1}z\|_{0,{\Delta_{T}}}}.

As for the cell contribution to the a posteriori estimate, we can retain, for each facet, the more favorable estimator depending on whether the facet $F$ belongs to an element in ${\mathcal{T}}_{h}^{1}$ or in ${\mathcal{T}}_{h}^{2}$ . By similar argument to the ones used for the element residual term, we have

(23)

-\sum_{F\in{{\mathcal{F}_{h}^{i}}}}\int_{F}\Lbrack a\partial_{\nu}u_{h}\Rbrack(z-\widehat{\Pi}_{h}z)\lesssim C(\Omega)\sqrt{\sum_{e\in{{\mathcal{F}_{h}^{i}}}}\varsigma_{F}^{2}|\mathbf{r}_{0}(F)|^{2}}.

The boundary terms are treated in the standard way for any $F\in{{\mathcal{F}_{h}^{b}}}$ and $F\subset\partial T$ ,

(24)

\begin{split}&\int_{F}(\lambda_{h}-a\partial_{\nu}u_{h})(z-\widehat{\Pi}_{h}z)\leq\|\lambda_{h}-a\partial_{\nu}u_{h}\|_{0,F}\|z-\widehat{\Pi}_{h}z\|_{0,F}\lesssim\|\lambda_{h}-a\partial_{\nu}u_{h}\|_{0,F}h_{F}^{1/2}|z|_{1,{\Delta_{T}}}.\end{split}

Therefore,

(25)

\sum_{F\in{{\mathcal{F}_{h}^{b}}}}\int_{F}(\lambda_{h}-a\partial_{\nu}u_{h})(z-\widehat{\Pi}_{h}z)\lesssim\left(\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}\|\lambda_{h}-a\partial_{\nu}u_{h}\|^{2}_{0,F}\right)^{1/2}\|z\|_{1,\Omega}\leq\left(\sum_{F\in{{\mathcal{F}_{h}^{b}}}}|\mathbf{r}_{1}(F)|^{2}\right)^{1/2}.

By Eq. 9, the last term can be bounded as

\int_{\Gamma}(g-u_{h})\zeta\leq\|g-u_{h}\|_{1/2,\Gamma}\|\zeta\|_{-1/2,\Gamma}\lesssim\|g-u_{h}\|_{1/2,\Gamma}.

Finally, since $g-u_{h}$ is orthogonal to $\widecheck{\Lambda}_{h}\subseteq\Lambda_{h}$ , we can use Lemma 3 of [15] to bound

(26)

\|g-u_{h}\|^{2}_{1/2,\Gamma}\lesssim\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}|g-u_{h}|_{1,F}^{2}=\sum_{F\in{{\mathcal{F}_{h}^{b}}}}|\mathbf{r}_{2}(F)|^{2}.

Combining all gives Eq. 18. This completes the proof of the theorem.

If $g\in H^{1}(\Gamma)$ , it is, therefore, natural, for the Lagrangian multiplier method, to define the following error indicator $\eta_{T}$ for each element $T\in{\mathcal{T}}_{h}$ , and estimator $\eta$ , by

(27)

\begin{split}\eta_{T}&=\sqrt{\varsigma_{T}^{2}|{\mathbb{r}}(T)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{i}}}\cap\partial T}\varsigma_{F}^{2}|\mathbf{r}_{0}(F)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{b}}}\cap\partial T}\left(|\mathbf{r}_{1}(F)|^{2}+|\mathbf{r}_{2}(F)|^{2}\right)},\\ \eta&=\sqrt{\sum_{T\in{\mathcal{T}}_{h}}\eta_{T}^{2}}.\end{split}

If $g$ is not in $H^{1}(\Gamma)$ , we can not get the full localization Eq. 26 of the residual $g-u_{h}$ on $\Gamma$ . We can however resort, for the two dimensional case, to [22, Theorem 2.2], and, for the three dimensional case, to [23, Lemma 3.1], which allow us to bound

\|g-u_{h}\|^{2}_{1/2,\Gamma}\lesssim\sum_{P\in\mathcal{N}^{b}_{h}}|g-u_{h}|^{2}_{1/2,{\Delta_{P}}}

where $\mathcal{N}^{b}_{h}$ is the set of nodes of the mesh ${\mathcal{T}}_{h}$ on $\Gamma$ and where for $P\in\mathcal{N}^{b}_{h}$ , ${\Delta_{P}}\subset\Gamma$ is the patch formed by the boundary facets sharing $P$ as a vertex. For those patches ${\Delta_{P}}$ for which $g|_{{\Delta_{P}}}\in H^{1}({\Delta_{P}})$ , $|g-u_{h}|_{1/2,{\Delta_{P}}}$ can be further bounded by $\underset{F\subset{\Delta_{P}}}{\sum}|\mathbf{r}_{2}(F)|^{2}$ . For the remaining patches, the $H^{1/2}({\Delta_{P}})$ semi-norm of the residual will have to be computed by evaluating the double integral involved in the definition of the fractional norm.

Remark 2.9.

Note that, in the implementation of the method, we do not explicitly use the splitting $\mathcal{T}_{h}=\mathcal{T}_{h}^{1}\cup\mathcal{T}_{h}^{2}$ , which is only needed for the theoretical analysis. Remark also that, for the elements adjacent to $\Gamma$ , for which $\rho^{-k}_{T}=\infty$ , we always have $\varsigma_{T}=C_{1}$ .

3 Application to stabilized methods for the imposition of boundary conditions

In engineering practice it is often advantageous to use a stabilized method instead of choosing the spaces so that the inf-sup condition is satisfied. In this section we show how the proposed framework can be adapted to two of the most well-known stabilized methods, namely the Barbosa-Hughes method [4] and the Nitsche’s method [37]. We assume for the sake of simplicity that $g\in H^{1}(\Gamma)$ . Both the final results and the arguments are in the same spirit as Theorem 2.7 above and therefore we only give sketches of the proofs.

3.1 Indicators for the Barbosa–Hughes method

The Barbosa–Hughes discrete problem reads: find $u_{h}\in V_{h}$ , $\lambda_{h}\in\Lambda_{h}$ such that for all $v_{h}\in V_{h}$ , $\mu_{h}\in\Lambda_{h}$ it holds that

(28)		$\displaystyle\int_{\Omega}a\nabla u_{h}\cdot\nabla v_{h}-\int_{\Gamma}\lambda_{h}v_{h}\pm\alpha\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}\int_{F}(a\partial_{\nu}u_{h}-\lambda_{h})(a\partial_{\nu}v_{h})=\int_{\Omega}fv_{h},$
(29)		$\displaystyle\int_{\Gamma}u_{h}\mu_{h}-\alpha\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}\int_{F}(a\partial_{\nu}u_{h}-\lambda_{h})\mu_{h}=\int_{\Gamma}g\mu_{h}.$

Here we use $\pm$ in front of the stabilization term in Eq. 28, to indicate that the analysis applies to both the symmetric and antisymmetric version of the method. The functional $L$ and $z,\zeta$ are defined as in the previous section. Similarly we have the following error representation by subtracting Eq. 28 and Eq. 29 from Eq. 1: for arbitrary $z_{h}\in V_{h}$ and $\zeta_{h}\in\Lambda_{h}$ it holds that

\begin{split}L(\delta)=&\int_{\Omega}a\nabla e\cdot\nabla z-\int_{\Gamma}\delta z+\int_{\Gamma}e\zeta\\ =&\int_{\Omega}a\nabla e\cdot\nabla(z-z_{h})-\int_{\Gamma}\delta(z-z_{h})+\int_{\Gamma}e(\zeta-\zeta_{h})+\int_{\Omega}a\nabla e\cdot\nabla z_{h}-\int_{\Gamma}\delta z_{h}+\int_{\Gamma}e\zeta_{h}\\ =&\int_{\Omega}a\nabla e\cdot\nabla(z-z_{h})-\int_{\Gamma}\delta(z-z_{h})+\int_{\Gamma}e(\zeta-\zeta_{h})-\alpha\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}\int_{F}(a\partial_{\nu}u_{h}-\lambda_{h})(\zeta_{h}\mp a\partial_{\nu}z_{h}).\end{split}

From Proposition 2.1, we have

(30)

\begin{split}L(\delta)=&\sum_{T\in{\mathcal{T}}_{h}}\int_{T}(f+{\nabla\cdot a\,\nabla}u_{h})(z-z_{h})-\sum_{F\in{{\mathcal{F}_{h}^{i}}}}\int_{F}\Lbrack a\partial_{\nu}u_{h}\Rbrack(z-z_{h})+\sum_{F\in{{\mathcal{F}_{h}^{b}}}}\int_{F}(\lambda_{h}-a\partial_{\nu}u_{h})(z-z_{h})\\ &+\int_{\Gamma}(g-u_{h})(\zeta-\zeta_{h})-\alpha\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}\int_{F}(a\partial_{\nu}u_{h}-\lambda_{h})(\zeta_{h}\mp a\partial_{\nu}z_{h}).\end{split}

We again set $z_{h}=\widehat{\Pi}_{h}z$ , $\zeta_{h}=0$ . The first three terms in Eq. 30 can be bounded using Eq. 20, Eq. 23 and Eq. 25. However, for the fourth term in Eq. 30, contrary to the previous case, we do not have that $u_{h}-g$ is orthogonal to the multiplier space, therefore Eq. 26 no longer holds. Instead we only have the following weaker bound [22, 23]. Recall that $\mathcal{N}^{b}_{h}$ denote the set of boundary vertices of the triangulation, and for each $P\in\mathcal{N}^{b}_{h}$ denote by ${\Delta_{P}}\subset\Gamma$ the patch formed by the boundary faces sharing $P$ as a vertex. We have

(31)

\|g-u_{h}\|_{1/2,\Gamma}^{2}\lesssim\sum_{P\in\mathcal{N}^{b}_{h}}|u_{h}-g|^{2}_{1/2,{\Delta_{P}}}+\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}^{-1}\|u_{h}-g\|^{2}_{0,F}.

We can further localize the term $|u_{h}-g|^{2}_{1/2,{\Delta_{P}}}$ . In order to do so we add and subtract $g^{P}_{h}\in\widecheck{V}_{h}|_{\Delta_{P}}$ , where $g^{P}_{h}$ is the $L^{2}({\Delta_{P}})$ projection onto the local space of continuous piecewise linears on the $(d-1)$ dimensional local mesh ${\mathcal{T}}_{h}|_{\Delta_{P}}$ , yielding

|g-u_{h}|^{2}_{1/2,{\Delta_{P}}}\lesssim|u_{h}-g^{P}_{h}|^{2}_{1/2,{\Delta_{P}}}+|g^{P}_{h}-g|^{2}_{1/2,{\Delta_{P}}},

which, combining with the inverse inequality, gives

|u_{h}-g^{P}_{h}|^{2}_{1/2,{\Delta_{P}}}\lesssim{h_{P}^{-1}}\|u_{h}-g^{P}_{h}\|^{2}_{0,{\Delta_{P}}}\simeq\sum_{F\subseteq{\Delta_{P}}}h_{F}^{-1}\|u_{h}-g^{P}_{h}\|^{2}_{0,F},

where $h_{P}=\max_{F\in{\Delta_{P}}}h_{F}$ (remark that the shape regularity of the mesh implies that for all $F\subseteq{\Delta_{P}}$ we have $h_{F}\simeq h_{P}$ ). Thanks to the fact that $g-g_{h}^{P}$ is orthogonal to the continuous piecewise linear functions, we have

|g^{P}_{h}-g|^{2}_{1/2,{\Delta_{P}}}\lesssim\sum_{F\subseteq{\Delta_{P}}}h_{F}|g^{P}_{h}-g|^{2}_{1,F}.

We also observe that, for $P$ a vertex of $F$

\|g-u_{h}\|^{2}_{0,F}\lesssim\|g-g^{P}_{h}\|^{2}_{0,F}+\|g^{P}_{h}-u_{h}\|^{2}_{0,F}\lesssim\sum_{F\subseteq{\Delta_{P}}}h_{F}|g^{P}_{h}-g|^{2}_{1,F}+\|u_{h}-g^{P}_{h}\|^{2}_{0,F}.

Combining these bounds we easily obtain

(32)

\begin{split}\int_{\Gamma}(g-u_{h})\zeta&\leq\|g-u_{h}\|_{1/2,\Gamma}\|\zeta\|_{-1/2,\Gamma}\lesssim\|g-u_{h}\|_{1/2,\Gamma}\\ &\lesssim\sqrt{\sum_{P\in\mathcal{N}^{b}_{h}}\sum_{F\subseteq{\Delta_{P}}}(h_{F}^{-1}\|u_{h}-g^{P}_{h}\|^{2}_{0,F}+h_{F}|g^{P}_{h}-g|_{1,F}^{2})}=\sqrt{\sum_{P\in\mathcal{N}^{b}_{h}}\sum_{F\subseteq{\Delta_{P}}}|\mathbf{r}(F,P)|^{2}},\end{split}

where, for $P$ a vertex of $F\subset{{\mathcal{F}_{h}^{b}}}$ we define

\mathbf{r}(F,P)=\sqrt{h_{F}^{-1}\|u_{h}-g^{P}_{h}\|^{2}_{0,F}+h_{F}|g^{P}_{h}-g|_{1,F}^{2}}\,.

Finally, we bound the additional term resulting from the stabilization, namely

(33)

\begin{split}\sum_{F\in{{\mathcal{F}}_{h}^{b}}}h_{F}\int_{F}(a\partial_{\nu}u_{h}-\lambda_{h})(a\partial_{\nu}(\widehat{\Pi}_{h}z))&\leq\sqrt{\sum_{F\in{{\mathcal{F}}_{h}^{b}}}h_{F}\|a\nabla u_{h}\cdot\nu-\lambda_{h}\|_{0,F}^{2}}\sqrt{\sum_{F\in{{\mathcal{F}}_{h}^{b}}}h_{F}\|a\partial_{\nu}(\widehat{\Pi}_{h}z)\|_{0,F}^{2}}\\ \lesssim&\sqrt{\sum_{F\in{{\mathcal{F}}_{h}^{b}}}h_{F}\|a\partial_{\nu}u_{h}-\lambda_{h}\|_{0,F}^{2}}.\end{split}

The last bound derives from a standard trace inequality on the element $T$ associated to the boundary face $F$ followed by an inverse inequality and an $H^{1}$ stability bound for $\widehat{\Pi}_{h}$

(34)

\begin{split}\|a\partial_{\nu}(\widehat{\Pi}_{h}z)\|_{0,F}&\lesssim\|\nabla(\widehat{\Pi}_{h}z)\|_{0,F}\lesssim h_{T}^{-1/2}\|\nabla(\widehat{\Pi}_{h}z)\|_{0,T}+h_{T}^{1/2}|\nabla(\widehat{\Pi}_{h}z)|_{1,T}\\ &\lesssim h_{T}^{-1/2}\|\nabla(\widehat{\Pi}_{h}z)\|_{0,T}\lesssim h_{T}^{-1/2}|z|_{1,{\Delta_{T}}},\end{split}

which, together with (9), yields

\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}\|a\partial_{\nu}(\widehat{\Pi}_{h}z)\|_{0,F}^{2}\lesssim\|z\|^{2}_{1,\Omega}\lesssim 1.

We then have

(35)

\alpha\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}\int_{F}(a\nabla u_{h}\cdot\nu-\lambda_{h})(a\partial_{\nu}(\widehat{\Pi}_{h}z))\lesssim\alpha\sqrt{\sum_{F\in{{\mathcal{F}_{h}^{b}}}}|\mathbf{r}_{1}(F)|^{2}}

where, we recall, $\mathbf{r}_{1}(F)=h^{1/2}_{F}\|a\partial_{\nu}u_{h}-\lambda_{h}\|_{0,F}$ .

Collecting the above bounds we obtain the a posteriori error estimate for the Barbosa–Hughes formulation Eq. 28–Eq. 29:

(36)

\|\lambda-\lambda_{h}\|^{2}_{-1/2,\Gamma}\lesssim\sum_{T\in{\mathcal{T}}_{h}}\varsigma_{T}^{2}|{\mathbb{r}}(T)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{i}}}}\varsigma_{F}^{2}|\mathbf{r}_{0}(F)|^{2}+(1+\alpha^{2})\sum_{F\in{{\mathcal{F}_{h}^{b}}}}|\mathbf{r}_{1}(F)|^{2}+\sum_{P\in\mathcal{N}^{b}_{h}}\sum_{F\subseteq{\Delta_{P}}}|\mathbf{r}(F,P)|^{2}.

3.2 Indicators for Nitsche’s method

Let us now consider Nitsche’s method, which reads: find $u_{h}\in V_{h}$ such that for all $v_{h}\in V_{h}$ , there holds

(37)

\begin{split}&\int_{\Omega}a\nabla u_{h}\cdot\nabla v_{h}-\int_{\Gamma}v_{h}(a\partial_{\nu}u_{h})\pm\int_{\Gamma}u_{h}(a\partial_{\nu}v_{h})+\gamma\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}^{-1}\int_{F}u_{h}v_{h}\\ =&\int_{\Omega}fv_{h}\pm\int_{\Gamma}g(a\partial_{\nu}v_{h})+\gamma\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}^{-1}\int_{F}gv_{h}.\end{split}

Following the work of Stenberg [41], which focuses on the Poisson equation but which is easily adapted to Eq. 1, Nitsche’s method is equivalent to a Barbosa-Hughes method with the choice $\Lambda_{h}=L^{2}(\Gamma)$ . The solution $u_{h},\lambda_{h}$ of Eq. 28-Eq. 29 with $\Lambda_{h}=L^{2}(\Gamma)$ verifies that $u_{h}$ solves Eq. 37 with $\gamma=\alpha^{-1}$ , and we have that, on $e\subset\Gamma$ ,

(38)

\lambda_{h}=a\partial_{\nu}u_{h}+\gamma h_{F}^{-1}(g-u_{h}).

Due to the equivalence, $L(\delta)$ has the same representation Eq. 30 with $\alpha$ replaced by $\gamma^{-1}$ . With the same choice of $z_{h}$ and $\zeta_{h}$ , the first and second terms can be bounded using Eq. 20 and Eq. 23 respectively. The fourth term $\int_{\Gamma}(g-u_{h})\zeta$ can be bounded using Eq. 32. For the remaining terms, observing that

\mathbf{r}_{1}(F)=h_{F}^{1/2}\|\lambda_{h}-a\partial_{\nu}u_{h}\|_{0,F}=h_{F}^{1/2}\|\gamma h_{F}^{-1}(g-u_{h})\|_{0,F}:=\gamma\mathbf{r}_{3}(F),

with

\mathbf{r}_{3}(F)=h_{F}^{-1/2}\|u_{h}-g\|_{0,F},

which, combining with Eq. 25 and Eq. 35, yields

(39)

\sum_{F\in{{\mathcal{F}_{h}^{b}}}}\int_{F}(\lambda_{h}-a\partial_{\nu}u_{h})(z-\widehat{\Pi}_{h}z)\mp\gamma^{-1}\sum_{F\in{{\mathcal{F}_{h}^{b}}}}h_{F}\int_{F}(a\partial_{\nu}u_{h}-\lambda_{h})(a\nabla(\widehat{\Pi}_{h}z)\cdot\nu)\lesssim\sqrt{\sum_{F\in{{\mathcal{F}_{h}^{b}}}}(1+\gamma^{2})|\mathbf{r}_{3}(F)|^{2}}.

Collecting all, we obtain the following a posteriori error bound for the normal flux computed using Nitsche’s method.

(40)

\|a\partial_{\nu}u-\lambda_{h}\|_{-1/2,\Gamma}\lesssim\sqrt{\varsigma_{T}^{2}\sum_{T\in{\mathcal{T}}_{h}}|{\mathbb{r}}(T)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{i}}}}\varsigma_{F}^{2}|\mathbf{r}_{0}(F)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{b}}}}(1+\gamma^{2})|\mathbf{r}_{3}(F)|^{2}+\sum_{P\in\mathcal{N}^{b}_{h}}\sum_{F\subseteq{\Delta_{P}}}|\mathbf{r}(F,P)|^{2}},

where $\lambda_{h}$ is given in Eq. 38.

Similarly as in Eq. 27, we can define the corresponding error indicator $\eta_{T}$ and $\eta$ for the Barbosa–Hughes and Nitsche’s methods.

4 Numerical experiments

In this section we demonstrate the performance of the proposed error estimator on some simple, yet significant, two dimensional test cases. Firstly, we demonstrate the action of the weight $\varsigma_{T}$ defined in Eq. 19 in the adaptive mesh refinement procedure independently of any particular problem. For simplicity, we fix $C_{1}=1$ . In the computation, we approximate $\rho_{T}$ by the following:

\rho_{T}\approx\min_{x\in\mathcal{N}_{\triangle_{T}}}d_{\Gamma}(x)

where $\mathcal{N}_{\triangle_{T}}$ is the set of all vertices on $\triangle_{T}$ .

We start with a $4$ by $4$ initial triangular mesh on a unit square domain, see Fig. 1(a). A total number of $7$ refinement steps are performed and the marking strategy identifies an element $K\in\mathcal{T}_{h}$ to be refined if

\varsigma_{T}>0.5\varsigma_{T,max},\quad\mbox{where}\,\varsigma_{T,max}=\max_{T\in\mathcal{T}_{h}}\varsigma_{T}.

Fig. 1 shows the meshes at various steps with $k=2$ and $C_{2}=1.0$ . It is easy to observe that significantly more refinements are placed near the boundary. Further experiments also show that the refinements on the boundary become more dominant if we decrease the value of $C_{2}$ or increase the order of $k$ .

Remark 4.1.

We note that, while the precise value of the best constants $C_{1}$ and $C_{2}$ in Lemma 2.5 is not known, it is possible to give an estimate of the ratio $C_{2}/C_{1}$ , in terms of the Poincaré constant for the patch ${\Delta_{T}}$ . Indeed, as $\widehat{\Pi}_{h}$ preserves polynomials of degree not greater than $k$ , we can write

\|u-\widehat{\Pi}_{h}u\|_{0,T}+h_{T}\|u-\widehat{\Pi}_{h}u\|_{1,T}\leq C_{1}h_{T}\inf_{p\in\mathbb{P}_{k}}|u-p|_{1,{\Delta_{T}}},

and then we can choose $C_{2}$ such that $C_{2}/C_{1}$ is the smallest constant $C_{k,\text{best}}$ for which

\inf_{p\in\mathbb{P}_{k}}|u-p|_{1,{\Delta_{T}}}\leq C_{k,\text{best}}h_{T}^{k}|u|_{k+1,{\Delta_{T}}}.

Such a constant may be estimated by recursively applying some upper bound for the Poincaré constant for the patch ${\Delta_{T}}$ , which can be obtained, for instance, by the approach of [42]. Its dependence on the polynomial degree $k$ can also be taken into account. For this choice to be the most effective, we would however need the upper bounds for $C_{k,\text{best}}$ to be sharp. If this is not the case, we observe that the true error might present some more or less pronounced oscillations. In our numerical tests, we tried several different values of $C_{2}$ . In all the cases considered, setting $C_{2}$ between $0.1$ to $1$ turns out to be a reasonable choice. See also Remark 4.1.

Refer to caption — Figure 1: Adaptive meshes based on $\varsigma_{T}$ with $k=2$ and $C_{2}=1.0$

4.1 Computation of the true error

In this subsection, we present two methods to compute the true error, i.e., $\|\lambda-\lambda_{h}\|_{-1/2,\Gamma}$ , for the purpose of comparison. From Eq. 11 and Eq. 12,

\|\lambda-\lambda_{h}\|_{-1/2,\Gamma}^{2}=|\nabla w|_{\Omega}^{2}\quad(\mbox{or }\left<\lambda-\lambda_{h},w\right>_{\Gamma}),

where $w\in H^{1}(\Omega)$ satisfies the following variational problem:

(41)

\int_{\Omega}\nabla w\cdot\nabla v=\int_{\Gamma}(\lambda-\lambda_{h})v,\quad\mbox{and}\quad\int_{\Gamma}w=0\quad\forall\,v\in H^{1}(\Omega).

Note that (41) is a pure Neumann problem. The compatibility of the solution is guaranteed since $\int_{\Gamma}(\lambda-\lambda_{h})=0$ for all aforementioned numerical methods. We approximate the true error in each refinement step using a two order higher finite element method on a finer mesh (compared to the mesh used in the adaptive procedure). We let $w_{h}\in V_{h}^{k+2}$ denote the Galerkin projection of $w$ on $V_{h}^{k+2}=\{v\in H^{1}(\Omega):v|_{T}\in P^{k+2}(T)\quad\forall\,T\in\tilde{\mathcal{T}}_{h}\}$ . Here $\tilde{\mathcal{T}}_{h}$ is the finer mesh. We then approximate the error by

(42)

\|\lambda-\lambda_{h}\|_{-1/2,\Gamma}^{2}\approx|\nabla w_{h}|_{\Omega}^{2}\quad(\text{or }\int_{\Gamma}(\lambda-\lambda_{h})w_{h}).

When $\lambda$ does not have enough regularity, using (42) to accurately compute the true error becomes infeasible as a very fine mesh is required to guarantee the accuracy. We therefore introduce another method to compute the true error by exploring properties of the wavelet decomposition. Indeed, it is known that, by expanding a function in $H^{-1/2}(\Gamma)$ based on a suitable wavelet basis, an equivalent $H^{-1/2}(\Gamma)$ norm can be computed by taking a weighted $L^{2}$ norm of the coefficient vector. The latter can be efficiently computed by applying a wavelet transform [18]. This only requires computations on $\partial\Omega$ , therefore we are able to compute the true error to a satisfactory accuracy even for low regularity $\lambda$ .

More precisely, given $v\in H^{-1/2}(\Gamma)$ , we aim at computing $\|v\|_{-1/2,\Gamma}$ . In order to do so, we consider the sequence of spaces $\{V_{j}\}_{j=0}^{\infty}$ such that $V_{j}\subset L^{2}(\Gamma)$ is the space of piecewise constant functions on the embedded uniform grid on $\Gamma$ with mesh size $|\Gamma|2^{-j}$ . We denote by $\{x^{j}_{k}\}_{k=0}^{2^{j}-1}$ , the nodes of the corresponding mesh, which we assume to be ordered counter-clock wise. For $v\in V_{j}$ , we can compute the vector $\mathbf{v}_{j}$ of length $2^{j}$

\mathbf{v}_{j}:=\{v_{jk}\}_{k=0}^{2^{j}-1}\quad\mbox{and}\quad v_{jk}=\dfrac{2^{j/2}}{|\Gamma|}\int_{x^{j}_{k}}^{x^{j}_{k+1}}v.

$\{v_{jk}\}_{k=0}^{2^{j}-1}$ is regarded as the coefficients of the $L^{2}(\Gamma)$ orthonormal bases consisting of the normalized characteristic functions on the elements of the grid. As $V_{j}\subset V_{j+1}$ , for all level $j$ we can decompose $v_{j+1}\in V_{j+1}$ as $v_{j+1}=v_{j}+d_{j}$ , with $v_{j}\in V_{j}$ obtained by applying a suitable oblique projector $P_{j}$ to $v_{j+1}$ . This gives us a telescopic expansion of all function in $V_{M}$ as $v_{M}=v_{0}+\sum_{j=0}^{M-1}d_{j}$ , and, passing to the limit as $M$ goes to infinity, of all functions in $L^{2}(\Gamma)$ as $v=v_{0}+\sum_{j=0}^{\infty}d_{j}$ . Given $\mathbf{v}_{j+1}$ , we can compute $\mathbf{v}_{j}:=\{v_{jk}\}_{k=0}^{2^{j}-1}$ and $\mathbf{d}_{j}:=\{d_{jk}\}_{k=0}^{2^{j}-1}$ (this last one being the vector of coefficients of $d_{j}$ with respect to a suitable basis for the space $W_{j}=(1-P_{j})V_{j+1}$ ), by applying a low-pass filter $h$ (strictly related with the projector $P_{j}$ ), and the band-pass filter $g=[1,-1]$ :

v_{jk}=\sum_{l=0}^{L}\dfrac{\sqrt{2}}{2}h(l)\,v_{j+1,2k+l}\quad\mbox{and}\quad d_{jk}=\sum_{l=0}^{1}\dfrac{\sqrt{2}}{2}g(l)\,v_{j+1,2k+l}=\dfrac{\sqrt{2}}{2}\left(v_{j+1,2k}-v_{j+1,2k+1}\right),

where $L+1$ is the length of the low-pass filter $h$ . In the above computation the function $v$ is considered as periodic, so that, when the index $2k+l>2^{j+1}-1$ , we extend the vector $\mathbf{v}_{j+1}$ as $v_{j+1,2^{j+1}+k}=v_{j+1,k},k\geq 0$ . For suitable choices of the low pass filter $h$ , the following norm equivalence holds for all $v\in H^{-1/2}(\Omega)$ ([19])

\|v\|_{-1/2,\Gamma}^{2}\simeq\|\mathbf{v}_{0}\|_{2}^{2}+\sum_{j=0}^{\infty}2^{-j}\|\mathbf{d}_{j}\|_{2}^{2},

where $\|\cdot\|_{2}$ denotes the Euclidean norm. In our experiments we choose the so called (2,2)-biorthogonal wavelet (see [18]), for which the low pass filter $h$ is

h=\dfrac{\sqrt{2}}{2}[3/128,-3/128,-11/64,11/64,1,1,11/64,-11/64,-3/128,3/128].

By choosing $M$ big enough and projecting $v$ onto $V_{M}$ (in our tests we use the $L^{2}$ orthogonal projection), we approximate the norm by

(43)

\|v\|_{-1/2,\partial\Omega}^{2}\approx\|\mathbf{v}_{0}\|_{2}^{2}+\sum_{j=0}^{M-1}2^{-j}\|\mathbf{d}_{j}\|_{2}^{2}.

4.2 Test results

Before presenting the results of our numerical tests, let us recall what the dependence of the error on the number of degrees of freedom is expected to be for an order $k$ method on either a uniform or a boundary concentrated mesh: letting $h$ denote the mesh size on the boundary and $N$ the total number of degrees of freedom, we have $h\simeq N^{-1/2}$ for uniform meshes, and $h\simeq N^{-1}|\log(N)|$ for boundary concentrated meshes. For a smooth solution, the error on the normal flux for optimal order $k$ method will behave like $h^{k}$ , that is, $N^{-k/2}$ for uniform grids and $N^{-k}|\log(N)|^{k}$ for boundary concentrated meshes.

To assess the performance of our estimator, we test it on the Lagrangian method without stabilization and on Nitsche’s method (the Barbosa-Hughes method being equivalent to the latter). Nitsche’s method with polynomial degree $k$ is optimal, i.e., it yields an order $k$ rate of convergence, on uniform meshes (see Table 1). For the Lagrangian method, the rate of convergence depends on the choice of the multiplier. We test two choices: discontinuous piecewise polynomials of order $k^{\prime}=k-2$ and continuous polynomials of order $k^{\prime}=k$ . Both choices yield inf-sup stable discretizations, yet they are both suboptimal (see Table 2): the first choice only provides, for the normal flux, an approximation of order at most $k-1/2$ , at the cost of using an order $k$ method in the bulk, while, in the presence of corners, the second only allows for an order 1 approximation of the normal flux, independently of $k$ , as it involves approximating a discontinuous function (the normal flux, in the presence of corners) by means of continuous functions. We point out that we are in no way advocating such choices as recommended methods for solving the problem considered (other choices for the multiplier, see Remark 1.1, allowing for optimality, are of course to be preferred for the actual computation of the flux). However, considering such suboptimal cases allows us to put the robustness of our method to the test, and to show that the refinement driven by our estimator can somehow make up for the lack of optimality.

Example 4.2.

In this example, we consider the Poisson equation on the unit square domain with right hand side and boundary data chosen so that the solution is the Franke function [25]

\begin{split}u(x,y)=&0.75\exp{\left(-(9x-2)^{2}/4-(9y-2)^{2}/4\right)}+0.75\exp{(-(9x+1)^{2}/49-(9y+1)/10)}\\ &+0.5\exp{(-(9x-7)^{2}/4-(9y-3)^{2}/4)}-0.2\exp{(-(9x-4)^{2}-(9y-7)^{2})}.\end{split}

This function has two peaks at $(2/9,2/9)$ and $(7/9,1/3)$ and one sink at $(4/9,7/9)$ .

We firstly test the convergence rate of the true error $\|\lambda-\lambda_{h}\|_{-1/2,\Gamma}$ on uniform meshes. The true error is computed using the aforementioned two methods. We denote by $E_{1}$ the error computed by (42) and by $E_{2}$ the error computed using the wavelet in Eq. 43 with $M=20$ . The problem (42) is solved on a finer uniform mesh with mesh size $h=1/64$ . Tables 1 and 2 show the convergence rates for $E_{1}$ . Observe that these are in agreement with the expected convergence rates given by the standard error estimates for the two methods, that is order $1$ (resp $2$ ) for Nitsche’s method with $k=1$ (resp. $k=2$ ), and order $3/2$ (resp. $1$ ) for the Lagrangian multiplier method with $k=2$ , $k^{\prime}=0$ , (resp. $k=2$ , $k^{\prime}=2$ ). From Tables 1 and Table 2, we also observe that the ratio between $E_{2}$ and $E_{1}$ is relatively stable (the fluctuation of the ratio is likely caused by the inaccurate computation of $E_{1}$ ). In particular, for Nitsche’s method, the ration $E_{2}/E_{1}$ remains close to $0.25$ for both orders. These results, therefore, confirm that $E_{2}$ is equivalent to the true error for both the Nitsche and Lagrangian multiplier methods.

Table 1: Example 4.2: Convergence rates for Nitsche’s method on uniform meshes

	$k=1$				$k=2$
h	$E_{1}$	rate	$E_{2}$	$E_{2}/E_{1}$	$E_{1}$	rate	$E_{2}$	$E_{2}/E_{1}$
1/8	3.35E-1	0.40	8.53E-2	0.25	2.86E-1	3.31	5.17E-2	0.18
1/16	1.73E-1	0.95	4.51E-2	0.26	3.19E-2	3.16	7.50E-3	0.23
1/32	8.66E-2	1.00	2.26E-2	0.25	4.69E-3	2.77	1.44E-3	0.30
1/64	4.33E-2	1.00	1.13E-2	0.26	2.51E-4	2.10	8.15E-5	0.32

Table 2: Example 4.2: Convergence rates for Lagrangian Multiplier method on uniform meshes

	$k=2,k^{\prime}=0$				$k=2,k^{\prime}=2$
h	$E_{1}$	rate	$E_{2}$	$E_{2}/E_{1}$	$E_{1}$	rate	$E_{2}$	$E_{2}/E_{1}$
1/8	4.58E-2	1.88	1.14E-2	0.25	9.83E-3	2.67	3.13E-2	3.18
1/16	1.43E-2	1.68	3.97E-3	0.28	2.65E-3	1.89	7.74E-3	2.92
1/32	4.80E-3	1.57	1.40E-3	0.29	1.18E-3	1.15	2.60E-3	2.19
1/64	5.62E-4	1.54	1.82E-4	0.32	5.88E-4	1.01	1.15E-3	1.96

We now test the adaptive mesh refinement (AMR) procedure for the Lagrangian method. In the adaptive procedure, we set the stopping criteria such that the total number of degree of freedoms (DOFs) less than $20,000$ . The marking strategy is set such that an element $T$ is marked to be refined if $\eta_{K}\geq 0.5\eta_{K,max}$ . In this example, we set $C_{2}=1.0$ . The initial mesh is set to be the $4\times 4$ mesh in Fig. 1(a). For comparison, we also perform the adaptive mesh refinemet procedure using the classical residual based error estimator (AMRc) without any dual weights . For the Lagrangian method, it is defined as

\eta_{classical}=\sqrt{\sum_{T\in{\mathcal{T}}_{h}}|{\mathbb{r}}(T)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{i}}}}|\mathbf{r}_{0}(F)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{b}}}}\left(|\mathbf{r}_{1}(F)|^{2}+h_{F}^{-1}\|g-u_{h}\|_{0,F}^{2}\right)}

and for the Nitsche’s method [17] is defined as

\eta_{classical}=\sqrt{\sum_{T\in{\mathcal{T}}_{h}}{\color[rgb]{0,0,1}\definecolor[named]{pgfstrokecolor}{rgb}{0,0,1}|}{\mathbb{r}}(T)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{i}}}}|\mathbf{r}_{0}(F)|^{2}+\sum_{F\in{{\mathcal{F}_{h}^{b}}}}\gamma^{2}h_{F}^{-1}\|g-u_{h}\|_{0,F}^{2}}.

It is well known that $\eta_{classical}$ is optimal in minimizing the energy norm of the error, i.e., $\|\nabla(u-u_{h})\|_{0,\Omega}$ . Note that comparing with Eq. 27, the $H^{1}$ norm in $\mathbf{r}_{2}(F)$ is reduced to $L^{2}$ norm on $\Gamma$ with adjusted weights for the Lagrangian method.

Fig. 2 shows the final meshes for the adaptive Lagrangian Multiplier method ( $k=2,k^{\prime}=0$ ) using respective $\eta_{classical}$ (left) and $\eta$ (right). It can be seen that the mesh generated by $\eta_{classical}$ has dense refinements around the interior peaks and sinks while the mesh generated by $\eta$ has more dense refinements near the boundary and almost completely ignore the peaks and sinks in the interior domain.

In the log-log plots Fig. 3a–Fig. 3c, we compare the convergence of true errors and estimators. The purpose of the convergence figures is to compare the two adaptive procedures using respectively the dual-weighted and the classical non-weighted error estimators. From Fig. 3a, we see that the error driven by $\eta$ converges faster than the one driven by $\eta_{classical}$ , which already has the order $N^{-1}$ with $N$ being the total number of DOFs. In comparison with rates attained by uniform refinement, that are provided in the Table 2 and Table 1, the relationship is that the rate obtained by $\eta_{classical}$ is higher or equal than the uniform approximation rate, and that the rate obtained by $\eta$ is higher than that obtained by $\eta_{classical}$ . More in detail, in Fig. 3a we display two reference straight lines: the slope $-1$ of the first line refers to the approximation rate in the energy norm that can be attained by the best approximation with order $k$ finite elements on a quasi uniform grid with $N$ degrees of freedom, which also provides an upper bound for the corresponding error of the normal flux. The slope of the second reference line is numerically evaluated by linear regression of the data set $(\log(N),\log(E))$ from the AMR with the proposed estimator $\eta$ . For this case its value is $\sim-1.5$ . For the figures thereafter, the same strategies will be used to present the reference slopes.

Fig. 3b shows that both methods display the same rate of convergence with respect to the number of DOFs on the boundary. However, for the same total number of DOFs, much more DOFs are located on the boundary by $\eta$ . More precisely, Fig. 3c shows that the ratio between the numbers of boundary DOFs and the total DOFs gradually gets higher for the meshes generated by $\eta$ in the AMR procedure.

We also test Example 4.2 using the Lagrangian Multiplier method with $k=2$ , $k^{\prime}=2$ and

(44)

\Lambda_{h}=\{\lambda\in C^{0}(\Gamma):\ u|_{{{F}}}\in\mathbb{P}_{2}({{F}}),\ \forall{{F}}\in{\mathcal{T}}_{h}|_{\Gamma}\}.

Since in this test the domain has corners, and, consequently, $\lambda$ is discontinuous, optimal approximation for the multiplier can not be achieved, as $\lambda_{h}\in C^{0}(\Gamma)$ . This also shows in Table 2 for the uniform refinement. In Fig. 4b, we observe that the mesh is densely refined around the corners which indicates that the error estimator $\eta$ successfully captures the error on the corners.

We now test the Nitsche’s method with $k=1$ and $k=2$ and set $\gamma=10$ in Eq. 37. Fig. 6 compares the final meshes generated using $\eta$ and $\eta_{classical}$ . We observe similar phenomena to that of the Lagrangian method, i.e., the mesh generated by $\eta_{classical}$ has dense refinement near the interior peaks and sinks while the mesh generated by $\eta$ has dense refinements almost all close to the boundary. The corresponding convergence rates of the true error and error estimators are plotted in Fig. 8. Again, for both orders, we observe significant improvements of the convergence rate comparing to the classical case.

For the Nitsche’s method of linear order, we also compare the performance with the boundary concentrated meshes proposed in [38] by Pfefferer and Winkler, which yield what is presently the best a priori error estimate for a non adaptive approximation of the normal flux. The boundary concentrated mesh has a fixed hierarchy structure, i.e., it has uniform mesh size $h^{2}$ on the boundary and $h\sqrt{\mbox{dist}(T,\Gamma)}$ for interior elements. We generate three such meshes in Fig. 7.

The corresponding log-log curve of $(N,E_{2})$ is plotted with legend $E_{2,PW}$ in Fig. 8 (see the top left figure). We observe that the error obtained with this mesh, which is adapted “a priori” to a good approximation of the boundary flux, is very close to the error obtained thanks to our error estimator $\eta$ . This is not surprising, as the solution of the problem is smooth.

To provide a more complete picture, in Fig. 9, we instead compare the performance of the error estimators $\eta_{classical}$ and $\eta$ , and of the related adaptive mesh refinements, in terms of the convergences in the energy error $\|\nabla(u-u_{h})\|_{0,\Omega}$ . The results confirm that $\eta_{classical}$ is optimal for the energy error, while, as it is to be expected, $\eta$ yields only a sub-optimal rate for the energy error in each of the tests.

Remark that, despite the fact that both versions of the Lagrangian method that we tested are, for different reasons, suboptimal with respect to the order $k$ of the bulk discretization, our tests show that the proposed error estimator allows to obtain a satisfactory approximation of the normal flux also for such methods.

For the remaining examples, to avoid a too large number of redundant tests, we then focus only on Nitsche’s method, which is, instead, optimal and which, we recall, is equivalent to the Barbosa-Hughes method. Moreover, we observe that he results displayed before in Table 1 and Fig. 8 for Nitsche’s method both confirm that $E_{2}$ can, after rescaling, serve as a good alternative to the more expensive $E_{1}$ in evaluating the true error. As for Nitsche’s method with $k=1$ , the ratio $E_{2}/E_{1}$ is stable around $0.25$ , in the remaining examples, we will use $E=4E_{2}$ as the true error.

Example 4.3.

In this example, we test a diffusion problem with variable diffusion coefficient. The diffusion coefficient is defined as $a=1.0+\sin^{2}\left(\pi\sqrt{x^{2}+y^{2}}\right)$ . And the functions $g$ and $f$ are defined such that the true solution $u$ has the following representation:

u(x,y)=\exp(-\alpha_{p}((x-x_{p})^{2}+(y-y_{p})^{2}))\mbox{ with }\alpha_{p}=200,x_{p}=0.2,y_{p}=0.2.

Note that this function has a strong peak at the point $(x_{p},y_{p})$ .

In the adaptive procedure, the stopping criteria is again set such that the total number of DOFs is less than $20,000$ . We test the Nitsche’s method for both the first and second orders with $C_{2}=1.0$ . For Example 4.3 with variable coefficient, we observe similar numerical behavior as in Example 4.2, see Fig. 10–Fig. 11. From the left two sub-figures of Fig. 11, we observe that in both cases the convergence rates for the true error using $\eta$ is almost double than that using $\eta_{classical}$ . In the example, our adaptive algorithm slightly outperforms the PW method. In the case $k=2$ we however observe visible oscillations for the true error. This is not in contrast with the theory. Indeed, the Galerkin method minimizes a discrete energy norm of the error which controls the error on the normal flux only up to a constant. Therefore, refining the mesh does not automatically yield a reduction in the error on the normal flux, particularly if measured, as in our case, in a norm that does not depend on the diffusion coefficient $a$ .

Example 4.4.

In this example, we test the L-shaped domain Poisson problem ( $a=1$ ) with a corner singularity and with an addition interior peak. The true solution has the following representation in polar coordinates:

u(r,\vartheta)=r^{\alpha}\sin(\alpha\vartheta)+\exp(-\alpha_{p}((x-x_{p})^{2}+(y-y_{p})^{2}))\in H^{5/3}(\Omega)

where $\alpha=2/3$ , $(\alpha_{p},x_{p},y_{p})$ is the same as in Example 4.3, and the $\Omega$ is the L-shaped domain, i.e., $\Omega=[-1,1]^{2}\setminus(0,1)\times(-1,0)$ .

In this test, we set $C_{2}=1$ and $C_{2}=0.1$ for the first and second order Nitsche’s method, respectively. The convergence rate on uniform meshes is firstly verified in Table 3. We recall that, according to the standard a priori error estimates for uniformly refined grids, the error for both the Lagrangian multiplier and the Nitsche’s method behaves like $h^{5/3-1}=h^{2/3}=N^{-1/3}$ .

Table 3: Example 4.4: Convergence rates on uniform meshes

	Nitsche $k=1$		Nitsche $k=2$
h	$E_{2}$	rate	$E_{2}$	rate
1.76E-1	1.10E-2	0.86	6.88E-2	2.53
8.84E-2	1.28E-2	3.10	1.62E-2	2.08
4.42E-2	8.84E-3	0.54	7.93E-3	1.03
2.21E-2	6.04E-3	0.55	5.08E-3	0.64
1.10E-2	4.07E-3	0.56	3.36E-3	0.59
5.52E-3	2.72E-3	0.57	2.23E-3	0.59

The final meshes obtained for the Nitsche’s method are given in Fig. 12 and the corresponding convergence results are provided in Fig. 13. We note that for this problem, even with the low regularity caused by boundary singularity, in both cases the true error $E$ driven by $\eta$ still doubles the convergence rates with respect to those of $\eta_{classical}$ .

In this example, in the presence of the corner singularity on the boundary, the adaptive method based on our error estimator shows significantly better performance than the PW method, that has uniform refinement on the boundary.

Comparing the performance of the estimator in the three examples we see that the adaptive procedure based on the dual wighted residual performs always better than the one based on the classical error estimator. If the solution is smooth, the results obtained by the AMR based on the new estimator are, in terms of error vs number of degrees of freedom, as good as the ones obtained by using boundary concentrated meshes (of course, in this case, this last method is cheaper, as the mesh is designed a priori and the problem is solved only once). Our adaptive method is particularly advantageous when the solution presents singularities on or close to the boundary (which boundary concentrated meshes, based on a priori analysis, cannot tackle efficiently).

References

[1] M. Akira, A mixed finite element method for boundary flux computation, Computer methods in applied mechanics and engineering, 57 (1986), pp. 239–243, https://doi.org/10.1016/0045-7825(86)90016-2.
[2] T. Apel, J. Pfefferer, and A. Rösch, Finite element error estimates on the boundary with application to optimal control, Mathematics of Computation, 84 (2015), pp. 33–70, https://doi.org/10.1090/S0025-5718-2014-02862-7.
[3] T. Apel, J. Pfefferer, and M. Winkler, Local mesh refinement for the discretization of Neumann boundary control problems on polyhedra, Mathematical Methods in the Applied Sciences, 39 (2016), pp. 1206–1232, https://doi.org/10.1002/mma.3566.
[4] H. J. C. Barbosa and T. J. R. Hughes, The finite element method with Lagrange multipliers on the boundary: circumventing the Babuška-Brezzi condition, Comput. Methods Appl. Mech. Engrg., 85 (1991), pp. 109–128, https://doi.org/10.1016/0045-7825(91)90125-P.
[5] R. Becker, E. Estecahandy, and D. Trujillo, Weighted marking for goal-oriented adaptive finite element methods, SIAM J. Numer. Anal., 49 (2011), pp. 2451–2469, https://doi.org/10.1137/100794298.
[6] R. Becker, G. Gantner, M. Innerberger, and D. Praetorius, Goal-oriented adaptive finite element methods with optimal computational complexity, arXiv preprint, (2021), arXiv:2101.11407.
[7] R. Becker, M. Innerberger, and D. Praetorius, Optimal convergence rates for goal-oriented FEM with quadratic goal functional, arXiv e-prints, (2020), arXiv:2003.13270.
[8] R. Becker, H. Kapp, and R. Rannacher, Adaptive finite element methods for optimal control of partial differential equations: basic concept, SIAM J. Control Optim., 39 (2000), pp. 113–132, https://doi.org/10.1137/S0363012999351097.
[9] R. Becker and R. Rannacher, A feed-back approach to error control in finite element methods: basic analysis and examples, East-West J. Numer. Math., 4 (1996), pp. 237–264.
[10] R. Becker and R. Rannacher, Weighted a posteriori error control in FE methods, Citeseer, 1996.
[11] R. Becker and R. Rannacher, An optimal control approach to a posteriori error estimation in finite element methods, Acta Numer., 10 (2001), pp. 1–102, https://doi.org/10.1017/S0962492901000010.
[12] C. Bernardi, Y. Maday, and A. T. Patera, Domain decomposition by the mortar element method, in Asymptotic and numerical methods for partial differential equations with critical parameters, Springer, 1993, pp. 269–286, https://doi.org/10.1007/978-94-011-1810-1_17.
[13] S. Berrone, A. Bonito, R. Stevenson, and M. Verani, An optimal adaptive fictitious domain method, Math. Comp., 88 (2019), pp. 2101–2134, https://doi.org/10.1090/mcom/3414.
[14] S. Bertoluzza, Local boundary estimates for the Lagrange multiplier discretization of a Dirichlet boundary value problem with application to domain decomposition, Calcolo, 43 (2006), pp. 121–149, https://doi.org/10.1007/s10092-006-0115-7.
[15] S. Bertoluzza, Analysis of a mesh-dependent stabilization for the three fields domain decomposition method, Numerische Mathematik, 133 (2016), pp. 1–36, https://doi.org/10.1007/s00211-015-0742-5.
[16] F. Brezzi, L. P. Franca, D. Marini, and A. Russo, Stabilization Techniques for Domain Decomposition Methods with Non-matching grids, Citeseer, 1997.
[17] E. Burman, C. He, and M. G. Larson, A posteriori error estimates with boundary correction for a cut finite element method, IMA Journal of Numerical Analysis, (2020), https://doi.org/10.1093/imanum/draa085.
[18] A. Cohen, I. Daubechies, and J.-C. Feauveau, Biorthogonal bases of compactly supported wavelets, Communications on Pure and Applied Mathematics, 45 (1992), pp. 485–560, https://doi.org/10.1002/cpa.3160450502.
[19] W. Dahmen, Stability of multiscale transformations, Journal of Fourier Analysis and Applications, 2 (1996), pp. 341–361.
[20] B. H. Dennis and G. S. Dulikravich, Simultaneous determination of steady temperatures and heat fluxes on surfaces of three dimensional objects using FEM, ASME-Publications-HTD, 369 (2001), pp. 259–268.
[21] K. Eriksson and C. Johnson, Adaptive finite element methods for parabolic problems. I. A linear model problem, SIAM J. Numer. Anal., 28 (1991), pp. 43–77, https://doi.org/10.1137/0728003.
[22] B. Faermann, Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary elements methods. Part I. The two-dimensional case, IMA journal of numerical analysis, 20 (2000), pp. 203–234, https://doi.org/10.1093/IMANUM/20.2.203.
[23] B. Faermann, Localization of the Aronszajn-Slobodeckij norm and application to adaptive boundary element methods Part II. The three-dimensional case, Numerische Mathematik, 92 (2002), pp. 467–499.
[24] M. Feischl, D. Praetorius, and K. G. van der Zee, An abstract analysis of optimal goal-oriented adaptivity, SIAM J. Numer. Anal., 54 (2016), pp. 1423–1448, https://doi.org/10.1137/15M1021982.
[25] R. Franke, A critical comparison of some methods for interpolation of scattered data, tech. report, Navel Postgraduate School Monterey CA, 1979.
[26] D. Gilbarg and N. S. Trudinger, Elliptic partial differential equations of second order, Classics in Mathematics, Springer-Verlag, Berlin, 2001. Reprint of the 1998 edition.
[27] M. B. Giles, M. G. Larson, M. Levenstam, and E. Süli, Adaptive error control for finite element approximations of the lift and drag coefficients in viscous flow, 1997.
[28] V. Girault and R. Glowinski, Error analysis of a fictitious domain method applied to a Dirichlet problem, Japan J. Indust. Appl. Math., 12 (1995), pp. 487–514, https://doi.org/10.1007/BF03167240.
[29] P. M. Gresho, R. L. Lee, R. L. Sani, M. K. Maslanik, and B. E. Eaton, The consistent Galerkin FEM for computing derived boundary quantities in thermal and or fluids problems, International Journal for Numerical Methods in Fluids, 7 (1987), pp. 371–394, https://doi.org/10.1002/fld.1650070406.
[30] M. Holst and S. Pollock, Convergence of goal-oriented adaptive finite element methods for nonsymmetric problems, Numer. Methods Partial Differential Equations, 32 (2016), pp. 479–509, https://doi.org/10.1002/num.22002.
[31] M. Holst, S. Pollock, and Y. Zhu, Convergence of goal-oriented adaptive finite element methods for semilinear problems, Comput. Vis. Sci., 17 (2015), pp. 43–63, https://doi.org/10.1007/s00791-015-0243-1.
[32] T. Horger, J. M. Melenk, and B. Wohlmuth, On optimal L2- and surface flux convergence in FEM, Computing and Visualization in Science, 16 (2013), pp. 231–246, https://doi.org/10.1007/s00791-015-0237-z.
[33] M. Innerberger and D. Praetorius, Instance-optimal goal-oriented adaptivity, Comput. Methods Appl. Math., 21 (2021), pp. 109–126, https://doi.org/10.1515/cmam-2019-0115.
[34] B. Khoromskij and J. Melenk, Boundary concentrated finite element methods, SIAM J. Numer. Anal., 41 (2004), pp. 1–36.
[35] M. G. Larson and A. Massing, $L^{2}$ -error estimates for finite element approximations of boundary fluxes, arXiv e-prints, (2014), arXiv:1401.6994.
[36] J. Nerg and J. Partanen, A simplified FEM based calculation model for 3-D induction heating problems using surface impedance formulations, IEEE transactions on magnetics, 37 (2001), pp. 3719–3722, https://doi.org/10.1109/20.952698.
[37] J. Nitsche, Über ein variationsprinzip zur lösung von Dirichlet-problemen bei verwendung von teilräumen, die keinen randbedingungen unterworfen sind, Abhandlungen aus dem Mathematischen Seminar der Universität Hamburg, 36 (1971), pp. 9–15.
[38] J. Pfefferer and M. Winkler, Finite element error estimates for normal derivatives on boundary concentrated meshes, SIAM Journal on Numerical Analysis, 57 (2019), pp. 2043–2073, https://doi.org/10.1137/18M1181341.
[39] T. Richter and T. Wick, Variational localizations of the dual weighted residual estimator, Journal of Computational and Applied Mathematics, 279 (2015), pp. 192–208, https://doi.org/10.1016/j.cam.2014.11.008.
[40] L. R. Scott and S. Zhang, Finite element interpolation of nonsmooth functions satisfying boundary conditions, Math. Comp., 54 (1990), pp. 483–493, https://doi.org/10.2307/2008497.
[41] R. Stenberg, On some techniques for approximating boundary conditions in the finite element method, J. Comput. Appl. Math., 63 (1995), pp. 139–148, https://doi.org/10.1016/0377-0427(95)00057-7.
[42] A. Veeser and R. Verfürth, Poincaré constants for finite element stars, IMA Journal of Numerical Analysis, 32 (2011), pp. 30–47, https://doi.org/10.1093/imanum/drr011.


(a) step 1	(b) step 3

(c) step 5	(d) step 7


(a) $k=1$ by $\eta_{classical}$	(b) $k=1,C_{2}=1.0$ by $\eta$

(c) $k=2$ by $\eta_{classical}$	(d) $k=2,C_{2}=0.1$ by $\eta$


(a) Nitsche $k=1$	(b) Nitsche $k=2$

(c) Lagrangian $k=2,k^{\prime}=0$	(d) Lagrangian $k=2,k^{\prime}=2$


(a) $k=1,\eta_{classical}$	(b) $k=1,\eta$

(c) $k=2,\eta_{classical}$	(d) $k=2,\eta$


(a) $k=1$ , $\eta_{classical}$	(b) $k=1,C_{2}=1$ , $\eta$

(c) $k=2$ , $\eta_{classical}$	(d) $k=2,C_{2}=0.1$ , $\eta$

An a posteriori error estimate of the outer normal derivative using dual weights††thanks: Submitted to the editors of SIAM Journal of Numerical Analysis. \fundingEB and CH were funded by the EPSRC grant EP/P01576X/1.