Error Estimator Using Higher Order FEM for an Interface Problem

Maharavo Randrianarivony^1,2

¹Simulation Unit, Personal Simulation and Design, Sankt Augustin, Germany.
²Address: Pappelweg 7, Zimmer 21, Sankt Augustin 53757, Germany.

DOI: 10.4236/am.2017.812127                

Abstract

A higher order finite element method is considered to treat an interface problem. The polynomial degree is allowed to be arbitrary but it is fixed for the FEM-variational formulation. We propose an error estimator which turns out to be efficient and reliable. We demonstrate upper and lower bounds of the error estimator with respect to the exact accuracy. For the transmission problem, the coefficients for the internal and external domains can be highly dissimilar. One major difficulty is the characteristic of the estimator at the interface. The a-posteriori error estimates can be computed very efficiently element by element. To corroborate the theoretical analysis, we report on a few numerical results.

Keywords

Estimator, A-Posteriori, Interface, Transmission, Higher Order, FEM, Residual

1. Introduction

In most aspects of numerical simulation, it is desirable to provide an approxi- mation to an unknown. But it is even more reliable to supply an additional information quantifying how accurate the approximation is. In particular for FEM (Finite Element Method), that additional value is exactly the purpose of the a-posteriori error estimator which is described subsequently. We begin by motivating the transmission problem that is based on the PBE (Poisson-Boltz- mann Equation) for the interaction of solute and solvent media which are respectively denoted by ${\Omega }^{int}$ and ${\Omega }^{ext}$ . The surface $\Gamma$ represents the solute-solvent [1] [2] [3] interface which is the molecular surface in the realistic case. The solvent is represented by a continuous dielectric medium while the solute is located inside the cavity ${\Omega }^{int}$ . In the sequel, the whole solute-solvent domain is denoted by $\Omega \mathrm{<mo>\; :\; </mo>}={\Omega }^{int}\cup {\Omega }^{ext}$ . The nonlinear PBE [4] admits the next general expression

$-\nabla \cdot \left(\epsilon \left(x\right)\nabla u\left(x\right)\right)+{\displaystyle \underset{i=1}{\overset{M}{\sum }}}{n}_i{q}_i{\bar{\kappa }}^2\left(x\right){\text{e}}^{-\beta q_{i}u\left(x\right)}=\frac{4\text{π}{e}_C^2}{k_{B}T}{\displaystyle \underset{i=1}{\overset{N_m}{\sum }}}{z}_i\delta \left(x-x_i\right),\forall x\in \Omega$ (1)

The charge positions $x_i$ are located in the strict interior of the molecular surface $\Gamma$ as illustrated in Figure 1(a) and Figure 1(b). The quantities $k_B\mathrm{,}T\mathrm{,}{z}_i\mathrm{,}\beta \mathrm{,}{q}_i\mathrm{,}{n}_i\mathrm{,}{e}_C$ are constants related to physical parameters while the

(a)(b)(c)

Figure 1. Motivation: (a) Ionic solution; (b) Connolly surface: composition of trimmed toroidal and spherical surfaces; (c) Tetrahedral mesh of the solvent.

unknown function is the dimensionless electrostatic potential u. The coefficients $\epsilon \left(x\right)$ and $\kappa \left(x\right)$ are space-dependent functions related to the dielectric value and the modified Debye-Hückle parameter. Those coefficients might be discontinuous between ${\Omega }^{int}$ and ${\Omega }^{ext}$ but the solution u is required to be continuous everywhere. In the case $M=2$ and under certain assumption on the parameters, the exponential term becomes a hyperbolic sine. By the Taylor expansion, one has $sinh\left(t\right)=t+\left(t^3/\mathrm{3\; <mo>\; !\; </mo>}\right)+\left(t^5/\mathrm{5\; <mo>\; !\; </mo>}\right)+\cdots$ . The linear PBE considers only the first term of the above expansion to obtain $-\nabla \cdot \left(\epsilon \left(x\right)\nabla u\left(x\right)\right)+\mu \left(x\right)u\left(x\right)=f\left(x\right)$ for $x\in \Omega$ which will be the purpose of this paper. We will focus on the FEM treatment of the linearized equation for which we investigate a-posteriori error estimates.

Before presenting our approach, related works and our previous results are in order. Verfürth has compiled a comprehensive study [5] about APEE (a-posteriori error estimator) for which he mainly treats piecewise linear FEM. Many different a-posteriori error estimators have been proposed for the Stokes problem [6] [7] [8] [9] for isotropic grids. In the context of anisotropic meshes [10] [11] [12] , there are a variety of APEEs for Poisson and reaction-diffusion problems [13] [14] . An article [15] by Creusé, Kunert and Nicaise presents a survey on the residual based error estimator on anisotropic grids for the Stokes equation. An interesting APEE for two and three dimensions as well as an anisotropic adaptive mesh refinement are also detailed in [16] . Basically, a-posteriori error estimators permit to evaluate the finite element errors without knowing the exact solution. That feature makes it possible to dynamically identify regions where one should have further refinements if the error there is too large. Therefore, adaptive refinements are mainly based on the quality of a-posteriori error estimators. Our approach in this paper follows the same spirit as the works in [17] [18] . For the Spectral Element Method, we find in [18] an APEE for the hp-case. That is, the mesh size h is allowed to be refined in some regions while the polynomial degrees p are also variable on different elements of the mesh. The hp-case does not require that the polynomial degree or the mesh size are fixed. That has been generalized in [17] to treat hp-FEM [19] for the Poisson problem where corner singularities are allowed. We have presented in [20] an APEE based on hierarchical space enrichment on anisotropic FEM which is combined with adaptive refinements. Boundary Element Method (BEM) is very efficient [21] [22] [23] [24] [25] as far as the linear PBE is concerned because of the existence of a fundamental solution providing an explicit kernel for the integral equation formulation. In addition, BEM requires only a discretization of the surface $\Gamma$ instead of the entire volume $\Omega$ (see Figure 1(c)). When treating nonlinear PBE, a solver on the volume $\Omega$ appears to be unavoidable. This paper can be viewed as a preliminary work toward nonlinear PBE. We are still reluctant to completely focus on the nonlinear PDE because the equation in (1) presents a real challenge related to the exponential nonlinear term on the left hand side and the nonsmooth Dirac functions on the right hand side. The only treatment of nonlinear PBE using BEM which we are aware of is [26] that is admittedly a very good approach. By inspecting that paper in detail, we realized that a solver in the volume $\Omega$ is also needed to construct an artificial fundamental solution. An integral equation solved by BEM is then used by applying that artificial fundamental solution in order to form a kernel. That is, a treatment exclusively on the surface $\Gamma$ without recourse to a solver in $\Omega$ is so far not sufficient. Holst [27] [28] [29] [30] is one of the most prominent specialists of PBE using FEM. His work seems to be extensively based on piecewise linear variational formulation. The Finite Difference Method (FDM) is also widely used in PBE. The main reason does not seem to be attributed to its numerical efficiency but rather to code availability and to reference or comparison purpose (see Section 1.1.2 of [31] ). An important component of PBE simulation is the geometric information because exact solutions of PBE are only known for very few simple geometries. Implementing a program for generating an SES (Solvent Excluded Surface) from nuclei coordinates and the Van-der-Waals radii of the atoms is not straightforward because a lot of geometric tasks [32] [33] [34] come into play. It is a long process to start from the nuclei coordinates till obtaining geometric data for computations. We have achieved a geometric task to process nuclei information in order to generate data for BEM as well as a mesh generation [35] from FEM as illustrated in Figure 1(c). Furthermore, a real chemical simulation by using wavelet BEM is described in [25] for the quantum computation. A wavelet BEM simulation using domain decomposition techniques was described in [36] which treats the case of ASM (Additive Schwarz Method). It was utilized as an efficient preconditioner for the wavelet single layer potential which is badly conditioned. A wavelet BEM formulation for computing apparent surface charge is documented in [37] for an interface problem. A simulation for chemical quantum computation using FEM is documented in [38] where we used nanotubes immersed in polymer matrices.

We consider in this work a higher order FEM formulation to treat the interface problem. That is, the polynomial degree, which is arbitrary but fixed, is uniformly constant in the entire FEM mesh. We examine the a-posteriori error estimator locally within each element and each edge. There are several types of edges: the interface edges, the interior edges, the exterior edges and the boundary edges. The difficulty for an estimator with respect to an interface edge is the discontinuous coefficients on the incident elements as well as the flux jump at the interface. In this work, we are more interested in elaborating mathematical theory than in chemical simulation. The error estimator can be efficiently computed element by element. We consider smooth load functions as right hand side of the equation. In addition to the theoretical investigation, we contribute about the numerical influence of the parameters appearing in the a-posteriori error estimator. We need numerical tests because the dependence of all various constants with respect to the problem parameters is not established theoretically. The next discussion is structured as follows. In the following section, we start by formulating the problem accurately and we introduce some important definitions as well as the expression of the estimator. That is followed by the analysis of the a-posteriori error estimator in Section 1. We report on some interesting practical results in the last section.

2. Problem Setting and A-Posteriori Estimators

This section describes the problem formulation, the introduction of the higher order FEM as well as the explicit expressions of the a-posteriori error estimators. We recall also some important results related to polynomial inverse estimates. We consider the transmission PDE:

$-{\epsilon }^{int}\Delta u\left(x\right)+{\mu }^{int}u\left(x\right)=f\left(x\right)\text{for}x\in {\Omega }^{int}$ (2)

$-{\epsilon }^{ext}\Delta u\left(x\right)+{\mu }^{ext}u\left(x\right)=f\left(x\right)\text{for}x\in {\Omega }^{ext}$ (3)

$\underset{\begin{array}{l}x\to x_0\\ x\in {\Omega }^{int}\end{array}}{lim}u\left(x\right)=\underset{\begin{array}{l}x\to x_0\\ x\in {\Omega }^{ext}\end{array}}{lim}u\left(x\right)\text{for}{x}_0\in \Gamma$ (4)

${\epsilon }^{int}\underset{\begin{array}{l}x\to x_0\\ x\in {\Omega }^{int}\end{array}}{lim}\langle \nabla u\left(x\right)\mathrm{,}n\left(x_0\right)\rangle -{\epsilon }^{ext}\underset{\begin{array}{l}x\to x_0\\ x\in {\Omega }^{ext}\end{array}}{lim}\langle \nabla u\left(x\right)\mathrm{,}n\left(x_0\right)\rangle =Q\left(x_0\right)\text{for}{x}_0\in \Gamma$ (5)

$u\left(x\right)=0\text{for}x\in \partial \Omega$ (6)

where $n\left(x_0\right)$ designates the normal vector at $x_0\in \Gamma$ pointing outward of $\Gamma =\partial {\Omega }^{int}$ . The load function $f\mathrm{<mo>\; :\; </mo>}\Omega \to ℝ$ and the flux jump $Q\mathrm{<mo>\; :\; </mo>}\Gamma \to ℝ$ are given while the interface $\Gamma$ and the boundary $\partial \Omega$ are polygons in ${ℝ}^2$ . We consider a mesh ${\mathbb{M}}_h$ of the entire domain $\Omega ={\Omega }^{int}\cup {\Omega }^{ext}$ such that the restrictions of the mesh ${\mathbb{M}}_h$ in ${\Omega }^{int}$ and ${\Omega }^{ext}$ are respectively denoted by ${\mathbb{M}}_h^{int}$ and ${\mathbb{M}}_h^{ext}$ . The mesh ${\mathbb{M}}_h$ is composed of triangles admitting the next properties:

• The intersection of two different elements $T_i\mathrm{,}{T}_j\in {\mathbb{M}}_h$ is either empty or a common node or a complete edge,

• We have the coverings:

$\Omega ={\displaystyle \underset{T\in {\mathbb{M}}_h}{\cup }}T\mathrm{,}{\Omega }^{int}={\displaystyle \underset{T\in {\mathbb{M}}_h^{int}}{\cup }}T\mathrm{,}{\Omega }^{ext}={\displaystyle \underset{T\in {\mathbb{M}}_h^{ext}}{\cup }}T$ (7)

• Every node of the interface $\Gamma$ is also a node of ${\mathbb{M}}_h$ ,

• All edges of the interface $\Gamma$ are edges of the mesh ${\mathbb{M}}_h$ .

For a triangle $T\in {\mathbb{M}}_h$ , we denote

$h\left(T\right)\mathrm{<mo>\; :\; </mo>}=\text{diameter}\left(T\right)=sup\left\{{\Vert x-y\Vert }_{R^2}\mathrm{,}x\mathrm{,}y\in T\right\}$ (8)

$\rho \left(T\right)\mathrm{<mo>\; :\; </mo>}=\text{supremum}\text{of}\text{the}\text{diameters}\text{of}\text{all}\text{balls}\text{contained}\text{in}T$ (9)

$\sigma \left(T\right):=h\left(T\right)/\rho \left(T\right)=\text{aspect}\text{ratio}\text{of}T$ (10)

${\mathcal{N}}_0\left(T\right)\mathrm{<mo>\; :\; </mo>}=\text{set}\text{of}\text{elements}\text{of}{\mathbb{M}}_h\text{sharing}\text{a}\text{vertex}\text{with}T$ (11)

${\mathcal{N}}_i\left(T\right)\mathrm{<mo>\; :\; </mo>}=\text{set}\text{of}\text{elements}\text{of}{\mathbb{M}}_h\text{sharing}\text{a}\text{vertex}\text{with}{T}^{\prime }\in {\mathcal{N}}_{i-1}\left(T\right)$ (12)

We use $\mathcal{N}\left(T\right)$ to denote ${\mathcal{N}}_i\left(T\right)$ for sufficiently large i. We assume quasi-uniformity in the sense that there exists a constant ${\rho }_0>0$ such that

$\rho \left(T\right)\le {\rho }_0<\infty ,\forall T\in {\mathbb{M}}_h$ (13)

We define the following mutually disjoint subsets of edges

${\mathbb{E}}_h^{\Gamma }\mathrm{<mo>\; :\; </mo>}=\text{set}\text{of}\text{edges}\text{of}{\mathbb{M}}_h\text{on}\text{the}\text{interface}\Gamma$ (14)

${\mathbb{E}}_h^0\mathrm{<mo>\; :\; </mo>}=\text{set}\text{of}\text{edges}\text{of}{\mathbb{M}}_h\text{on}\text{the}\text{boundary}\partial \Omega \mathrm{,}$ (15)

${\mathbb{E}}_h^{int}\mathrm{<mo>\; :\; </mo>}=\text{set}\text{of}\text{edges}\text{of}{\mathbb{M}}_h^{int}\text{which}\text{are}\text{not}\text{included}\text{in}\Gamma \mathrm{,}$ (16)

${\mathbb{E}}_h^{ext}\mathrm{<mo>\; :\; </mo>}=\text{set}\text{of}\text{edges}\text{of}{\mathbb{M}}_h^{ext}\text{which}\text{are}\text{included}\text{neither}\text{in}\partial \Omega \text{nor}\text{in}\Gamma \mathrm{.}$ (17)

Note that an edge of ${\mathbb{E}}_h^{int}$ may have an endpoint in $\Gamma$ . Likewise, an edge of ${\mathbb{E}}_h^{ext}$ is allowed to have an endpoint in $\Gamma$ or $\partial \Omega$ . We introduce in addition the set of all edges

${\mathbb{E}}_h\mathrm{<mo>\; :\; </mo>}={\mathbb{E}}_h^{\Gamma }\cup {\mathbb{E}}_h^0\cup {\mathbb{E}}_h^{int}\cup {\mathbb{E}}_h^{ext}$ (18)

For an edge $e\in {\mathbb{E}}_h$ , we denote

$h\left(e\right)\mathrm{<mo>\; :\; </mo>}=\text{length}\left(e\right)=sup\left\{{\Vert x-y\Vert }_{R^2}\mathrm{,}x\mathrm{,}y\in e\right\}$ (19)

$\mathcal{N}\left(e\right)\mathrm{<mo>\; :\; </mo>}=\text{set}\text{of}\text{elements}\text{of}{\mathbb{M}}_h\text{having}e\text{as}\text{a}\text{side}$ (20)

$n\left(e\right)\mathrm{<mo>\; :\; </mo>}=\text{unit}\text{normal}\text{vector}\text{orthogonal}\text{to}e$ (21)

The direction of the normal vector $n\left(e\right)$ is outward $\Gamma$ if the edge $e\in {\mathbb{E}}_h^{\Gamma }$ while it is pointed toward the exterior of $\Omega$ if the edge $e\in {\mathbb{E}}_h^0$ . For all other edges in ${\mathbb{E}}_h^{int}\cup {\mathbb{E}}_h^{ext}$ , the normal vectors $n\left(e\right)$ are pointed in an arbitrary but fixed orientation.

For any triangle T, the affine invertible mapping from the unit reference

$\hat{T}\mathrm{<mo>\; :\; </mo>}=\left\{x=\left(x\mathrm{,}y\right)\in R^2\mathrm{<mo>\; :\; </mo>0}\le x\le \mathrm{1,}0\le y\le \mathrm{1,}0\le x+y\le 1\right\}$ (22)

onto T is denoted by $F_T\mathrm{<mo>\; :\; </mo>}\hat{T}\mathrm{<mo>\; :\; </mo>}\to T$ in which

$F_{T}x=Bx+b\text{where}det\left(B\right)=\mathcal{O}\left(h^2\left(T\right)\right)\mathrm{,}det\left(B^{-1}\right)=\mathcal{O}\left(h^{-2}\left(T\right)\right)$ (23)

That allows one to derive results on the unit reference triangle $\hat{T}$ and to carry them over to any element T in the original mesh ${\mathbb{M}}_h$ . We use the standard definitions of the Sobolev spaces for ${\mathbb{L}}^2\left(\Omega \right)$ , ${ℍ}^1\left(\Omega \right)$ and ${ℍ}_0^1\left(\Omega \right)$ . From here onward, we use the usual shorthand $X\precsim Y$ if there is a constant c such that $X\le cY$ in which c is independent on h and p. In addition, $X\simeq Y$ amounts to $X\precsim Y\precsim X$ .

We want to consider now the Galerkin variational formulation. Denote the restriction of the solution u to the interface problem in ${\Omega }^{int}$ and ${\Omega }^{ext}$ by $u^{int}$ and $u^{ext}$ respectively. Due to the Green identity we have

${\epsilon }^{int}{\displaystyle {\int }_{{\Omega }^{int}}}\nabla u^{int}\cdot \nabla v^{int}+{\mu }^{int}{\displaystyle {\int }_{{\Omega }^{int}}}{u}^{int}{v}^{int}={\displaystyle {\int }_{{\Omega }^{int}}}f{v}^{int}+{\epsilon }^{int}{\displaystyle {\int }_{\Gamma }}\frac{\partial u^{int}}{\partial n}{v}^{int}$ (24)

${\epsilon }^{ext}{\displaystyle {\int }_{{\Omega }^{ext}}}\nabla u^{ext}\cdot \nabla v^{ext}+{\mu }^{int}{\displaystyle {\int }_{{\Omega }^{ext}}}{u}^{ext}{v}^{ext}={\displaystyle {\int }_{{\Omega }^{ext}}}f{v}^{ext}-{\epsilon }^{ext}{\displaystyle {\int }_{\Gamma }}\frac{\partial u^{ext}}{\partial n}{v}^{ext}$ (25)

We will denote the piecewise constant function defined on $\Omega$ :

${\epsilon }^{\Omega }\left(x\right)\mathrm{<mo>\; :\; </mo>}=\{\begin{array}{l}{\epsilon }^{int}\text{if}x\in {\Omega }^{int}\hfill \\ {\epsilon }^{ext}\text{if}x\in {\Omega }^{ext}\hfill \end{array}{\mu }^{\Omega }\left(x\right)\mathrm{<mo>\; :\; </mo>}=\{\begin{array}{l}{\mu }^{int}\text{if}x\in {\Omega }^{int}\hfill \\ {\mu }^{ext}\text{if}x\in {\Omega }^{ext}\hfill \end{array}$ (26)

The sum of (24) and (25) for every $v\in {ℍ}_0^1\left(\Omega \right)$ yields

${\displaystyle {\int }_{\Omega }}{\epsilon }^{\Omega }\nabla u\cdot \nabla v+{\displaystyle {\int }_{\Omega }}{\mu }^{\Omega }uv={\displaystyle {\int }_{\Omega }}fv+{\displaystyle {\int }_{\Gamma }}\left({\epsilon }^{int}\frac{\partial u^{int}}{\partial n}-{\epsilon }^{ext}\frac{\partial u^{ext}}{\partial n}\right)v$ (27)

Taking into account the flux jump (5) in the transmission equation, we have

$Q={\epsilon }^{int}\langle n\mathrm{,}\nabla u^{int}\rangle -{\epsilon }^{ext}\langle n\mathrm{,}\nabla u^{ext}\rangle ={\epsilon }^{int}\frac{\partial u^{int}}{\partial n}-{\epsilon }^{ext}\frac{\partial u^{ext}}{\partial n}$ (28)

The Galerkin weak form is therefore

${\displaystyle {\int }_{\Omega }}{\epsilon }^{\Omega }\nabla u\cdot \nabla v+{\displaystyle {\int }_{\Omega }}{\mu }^{\Omega }uv={\displaystyle {\int }_{\Omega }}fv+{\displaystyle {\int }_{\Gamma }}Qv$ (29)

For a fixed polynomial degree $p\ge 1$ , the finite dimensional space is

${\mathbb{V}}_h\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}=\left\{v\in {\mathcal{C}}^0\left(\Omega \right)\mathrm{<mo>\; :\; </mo>}{v|}_T\in {\mathbb{P}}_p\mathrm{,}\forall T\in {\mathbb{M}}_h\right\}\text{in}\text{which}$ (30)

${\mathbb{P}}_p\mathrm{<mo>\; :\; </mo>}=\text{span}\left\{x^n{y}^m\mathrm{<mo>\; :\; </mo>0}\le n+m\le p\right\}$ (31)

The discrete Galerkin approximation is to search for $u_h\in {\mathbb{V}}_h\left(\Omega \right)$ such that

${\displaystyle {\int }_{\Omega }}{\epsilon }^{\Omega }\nabla u_h\cdot \nabla v_h+{\displaystyle {\int }_{\Omega }}{\mu }^{\Omega }{u}_h{v}_h={\displaystyle {\int }_{\Omega }}f{v}_h+{\displaystyle {\int }_{\Gamma }}Q{v}_h,\forall v_h\in {\mathbb{V}}_h\left(\Omega \right)$ (32)

Introduce the bilinear form

$\mathcal{A}\left(v\mathrm{,}w\right)\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_{\Omega }}{\epsilon }^{\Omega }\nabla v\cdot \nabla w+{\displaystyle {\int }_{\Omega }}{\mu }^{\Omega }vw$ (33)

In order to express the a-posteriori estimators, we assume that the appro- ximated solution $u_h$ on the current discretization ${\mathbb{M}}_h$ is available and we consider a parameter $\alpha \in \left[\mathrm{0,1}\right]$ . The 1D-weight on $\hat{e}=\left[0,1\right]$ and the 2D-weight on $\hat{T}$ are respectively

${\omega }_{\hat{e}}\left(t\right)=t\left(1-t\right)\text{and}{\omega }_{\hat{T}}\left(x\right)\simeq \text{distance}\left(x\mathrm{,}\partial \hat{T}\right)$ (34)

For a general edge e and triangle T, transformations from the reference elements are used to define ${\omega }_e$ and ${\omega }_T$ . For an interior element $T\in {�}_h^{int}$ , the estimator is defined as

${\left({\eta }_{\alpha \mathrm{,}T}^{int}\right)}^2\mathrm{<mo>\; :\; </mo>}=\frac{h{\left(T\right)}^2}{p^2}{\Vert \left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}^2\left(T\right)}^2$ (35)

where $f_T$ designates the ${\mathbb{L}}_2\left(T\right)$ -projection of the load function $f$ onto the element T. The expression for an exterior element $T\in {\mathbb{M}}_h^{ext}$ is similar:

${\left({\eta }_{\alpha \mathrm{,}T}^{ext}\right)}^2\mathrm{<mo>\; :\; </mo>}=\frac{h{\left(T\right)}^2}{p^2}{\Vert \left(f_T+{\epsilon }^{ext}\Delta u_h^{ext}-{\mu }^{ext}{u}_h^{ext}\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}^2\left(T\right)}^2$ (36)

For an interface edge $e\in {\mathbb{E}}_h^{\Gamma }$ , one introduces

${\left({\eta }_{\alpha \mathrm{,}e}^{\Gamma }\right)}^2\mathrm{<mo>\; :\; </mo>}=\frac{h\left(e\right)}{2p}{\Vert \left(Q_e-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n\left(e\right)}+{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}^2\left(e\right)}^2$ (37)

where $Q_e$ is the ${\mathbb{L}}_2\left(e\right)$ -projection of the flux jump Q onto the edge e. The estimator for an interior edge $e\in {\mathbb{E}}_h^{int}$ having a normal vector $n\left(e\right)$ is defined as

${\left({\eta }_{\alpha \mathrm{,}e}^{int}\right)}^2\mathrm{<mo>\; :\; </mo>}=\frac{h\left(e\right)}{2p}{\Vert 〚{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n\left(e\right)}〛{\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}^2\left(e\right)}^2$ (38)

where $〚t〛$ stands for the jump of t when evaluated from the two elements incident upon e. The orientation of the jump is irrelevant because one takes its square in the ${\mathbb{L}}^2$ -norm. The estimator corresponding to an exterior edge $e\in {\mathbb{E}}_h^{ext}$ is

${\left({\eta }_{\alpha \mathrm{,}e}^{ext}\right)}^2\mathrm{<mo>\; :\; </mo>}=\frac{h\left(e\right)}{2p}{\Vert 〚{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}〛{\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}^2\left(e\right)}^2$ (39)

Since one needs computable local estimators for an element-by-element computation, an interior element $T\in {\mathbb{M}}_h^{int}$ is introduced

${\left({\eta }_{\alpha \mathrm{,}T}^{loc}\right)}^2\mathrm{<mo>\; :\; </mo>}={\left({\eta }_{\alpha \mathrm{,}T}^{int}\right)}^2+{\displaystyle \underset{e\subset \partial T\mathrm{,}e\in {\mathbb{E}}_h^{int}}{\sum }}{\left({\eta }_{\alpha \mathrm{,}e}^{int}\right)}^2+{\displaystyle \underset{e\subset \partial T\mathrm{,}e\in {\mathbb{E}}_h^{\Gamma }}{\sum }}{\left({\eta }_{\alpha \mathrm{,}e}^{\Gamma }\right)}^2$ (40)

Likewise, for an exterior element $T\in {\mathbb{M}}_h^{ext}$ , the local estimator reads:

${\left({\eta }_{\alpha \mathrm{,}T}^{loc}\right)}^2\mathrm{<mo>\; :\; </mo>}={\left({\eta }_{\alpha \mathrm{,}T}^{ext}\right)}^2+{\displaystyle \underset{e\subset \partial T\mathrm{,}e\in {\mathbb{E}}_h^{ext}}{\sum }}{\left({\eta }_{\alpha \mathrm{,}e}^{ext}\right)}^2+{\displaystyle \underset{e\subset \partial T\mathrm{,}e\in {\mathbb{E}}_h^{\Gamma }}{\sum }}{\left({\eta }_{\alpha \mathrm{,}e}^{\Gamma }\right)}^2$ (41)

The local estimators add up to the global estimator:

${\eta }_{\alpha }^2\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{T\in {\mathbb{M}}_h}{\sum }}{\left({\eta }_{\alpha \mathrm{,}T}^{loc}\right)}^2$ (42)

One has the following polynomial inverse estimates and extension properties. The descriptions of the next lemmas are found in [17] [19] [39] [40] [41] [42] .

Lemma 1. Given $\alpha \mathrm{,}\beta$ such that $-1<\alpha <\beta$ and some $\delta \in \left[\mathrm{0,1}\right]$ . For every univariate polynomial ${\pi }_p$ of degree $p\ge 1$ on the 1D reference element $\hat{e}=\left[0,1\right]$ , one has

${\displaystyle {\int }_0^1}{\omega }_{\hat{e}}\left(t\right){\left[{{\pi }^{\prime }}_p\left(t\right)\right]}^2\text{d}t\le C_1{p}^2{\displaystyle {\int }_0^1}{\omega }_{\hat{e}}\left(t\right){\left[{\pi }_p\left(t\right)\right]}^2\text{d}t$ (43)

${\displaystyle {\int }_0^1}{\omega }_{\hat{e}}^{\alpha }\left(t\right){\left[{\pi }_p\left(t\right)\right]}^2\text{d}t\le C_2{p}^{2\left(\beta -\alpha \right)}{\displaystyle {\int }_0^1}{\omega }_{\hat{e}}^{\beta }\left(t\right){\left[{\pi }_p\left(t\right)\right]}^2\text{d}t$ (44)

${\displaystyle {\int }_0^1}{\omega }_{\hat{e}}^{2\delta }\left(t\right){\left[{{\pi }^{\prime }}_p\left(t\right)\right]}^2\text{d}t\le C_3{p}^{2\left(2-\delta \right)}{\displaystyle {\int }_0^1}{\omega }_{\hat{e}}^{\delta }\left(t\right){\left[{\pi }_p\left(t\right)\right]}^2\text{d}t$ (45)

On the 2D reference element (22), one has for every bivariate polynomial ${\pi }_p$ of degree $p\ge 1$

${\displaystyle {\int }_{\hat{T}}}{\omega }_{\hat{T}}\left(x\mathrm{,}y\right){\left|\nabla {\pi }_p\left(x\mathrm{,}y\right)\right|}^2\text{d}x\text{d}y\le C_1{p}^2{\displaystyle {\int }_{\hat{T}}}{\omega }_{\hat{T}}\left(x\mathrm{,}y\right){\left[{\pi }_p\left(x\mathrm{,}y\right)\right]}^2\text{d}x\text{d}y$ (46)

${\displaystyle {\int }_{\hat{T}}}{\omega }_{\hat{T}}^{\alpha }\left(x\mathrm{,}y\right){\left[{\pi }_p\left(x\mathrm{,}y\right)\right]}^2\text{d}x\text{d}y\le C_2{p}^{2\left(\beta -\alpha \right)}{\displaystyle {\int }_{\hat{T}}}{\omega }_{\hat{T}}^{\beta }\left(x\mathrm{,}y\right){\left[{\pi }_p\left(x\mathrm{,}y\right)\right]}^2\text{d}x\text{d}y$ (47)

${\displaystyle {\int }_{\hat{T}}}{\omega }_{\hat{T}}^{2\delta }\left(x\mathrm{,}y\right){\left|\nabla {\pi }_p\left(x\mathrm{,}y\right)\right|}^2\text{d}x\text{d}y\le C_3{p}^{2\left(2-\delta \right)}{\displaystyle {\int }_{\hat{T}}}{\omega }_{\hat{T}}^{\delta }\left(x\mathrm{,}y\right){\left[{\pi }_p\left(x\mathrm{,}y\right)\right]}^2\text{d}x\text{d}y$ (48)

The constants are $C_1=C_1\left(\alpha ,\beta \right)$ , $C_2=C_2\left(\alpha ,\beta \right)$ , $C_3=C_3\left(\delta \right)$ which do not depend on p.

Lemma 2. Consider a univariate polynomial $\pi$ of degree $p\ge 1$ defined on the unit reference interval $\hat{e}=\left[0,1\right]$ and a parameter $\mathrm{0\; <mo>\; <\; </mo>}\gamma \le 1$ . There exists a bivariate extension $v\in {ℍ}^1\left(\hat{T}\right)$ defined on the reference triangle $\hat{T}$ from (22) such that it has the next properties w.r.t. the weight ${\omega }_{\hat{e}}^{\alpha }$ :

${v|}_{\hat{e}}=\pi {\omega }_{\hat{e}}^{\alpha }\text{and}{v|}_{\left(\partial \hat{T}\backslash \hat{e}\right)}\equiv 0$ (49)

${\Vert v\Vert }_{{\mathbb{L}}_2\left(\hat{T}\right)}^2\le C\left(\alpha \right)\gamma {\Vert \pi {\omega }_{\hat{e}}^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(\hat{e}\right)}^2$ (50)

${\Vert \nabla v\Vert }_{{\mathbb{L}}_2\left(\hat{T}\right)}^2\le C\left(\alpha \right)\left[\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}\right]{\Vert \pi {\omega }_{\hat{e}}^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(\hat{e}\right)}^2$ (51)

3. Investigation of the Estimators

Theorem 1. Let $e$ be an interface edge in ${\mathbb{E}}_h^{\Gamma }$ . Define for the weight exponent $\alpha$ :

${\sigma }_e\mathrm{<mo>\; :\; </mo>}=\left(Q_e-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n\left(e\right)}+{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}\right){\omega }_e^{\alpha }$ (52)

We have for any $\gamma$ as in (50) and (51) the next estimate using the patch $\mathcal{N}\left(e\right)$ :

$\begin{array}{c}{\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}\le {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\{{\Vert f+{\epsilon }^T\Delta u_h-{\mu }^T{u}_h\Vert }_{{\mathbb{L}}_2\left(T\right)}\sqrt{\gamma h\left(T\right)}\\ +{\left|u-u_h\right|}_{{ℍ}^1\left(T\right)}\sqrt{\frac{1}{h\left(T\right)}\left(\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}\right)}\\ +{\Vert u-u_h\Vert }_{{\mathbb{L}}_2\left(T\right)}\sqrt{\gamma h\left(T\right)}\}+{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}\end{array}$ (53)

Under the assumption that $h^2\le Cp$ , we have the following bound:

$\begin{array}{c}{\eta }_{\alpha \mathrm{,}e}^{\Gamma }\le {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}{C}_{p\mathrm{,}\gamma }{\Vert u-u_h\Vert }_{{ℍ}^1\left(T\right)}+C_{\gamma }{\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\frac{h\left(T\right)}{\sqrt{p}}{\Vert f-f_T\Vert }_{{\mathbb{L}}_2\left(T\right)}\\ +\sqrt{\frac{h\left(e\right)}{p}}{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2\end{array}$ (54)

Proof. Let $T_e^{int}\in {\mathbb{M}}_h^{int}$ and $T_e^{ext}\in {\mathbb{M}}_h^{ext}$ be the two elements which are incident upon e. Designate by $\widetilde{{\sigma }_e}$ the extension as in (49) of the polynomial

${\pi }_p\mathrm{<mo>\; :\; </mo>}=Q_e-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n\left(e\right)}+{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}$ (55)

over the whole patch $\mathcal{N}\left(e\right)=T_e^{int}\cup T_e^{ext}$ such that we have the restriction ${\left(\widetilde{{\sigma }_e}\right)|}_e={\pi }_p{\omega }_e^{\alpha }={\sigma }_e$ and such that we have the boundary value ${\left(\widetilde{{\sigma }_e}\right)|}_{\partial \mathcal{N}\left(e\right)}=0$ . An application of the Green identity, which describes the partial integration of a function with respect to a domain and its boundary, on $T_e^{int}$ and $T_e^{ext}$ yields

$\begin{array}{l}{\displaystyle {\int }_{T_e^{int}}}\left(-{\epsilon }^{int}\Delta u^{int}+{\mu }^{int}{u}^{int}\right)\widetilde{{\sigma }_e}\\ ={\epsilon }^{int}{\displaystyle {\int }_{T_e^{int}}}\nabla u^{int}\cdot \nabla \widetilde{{\sigma }_e}-{\epsilon }^{int}{\displaystyle {\int }_e}\frac{\partial u_h^{int}}{\partial n\left(e\right)}+{\mu }^{int}{\displaystyle {\int }_{T_e^{int}}}{u}^{int}\widetilde{{\sigma }_e}\end{array}$ (56)

$\begin{array}{l}{\displaystyle {\int }_{T_e^{ext}}}\left(-{\epsilon }^{ext}\Delta u^{ext}+{\mu }^{ext}{u}^{ext}\right)\widetilde{{\sigma }_e}\\ ={\epsilon }^{ext}{\displaystyle {\int }_{T_e^{ext}}}\nabla u^{ext}\cdot \nabla \widetilde{{\sigma }_e}+{\epsilon }^{ext}{\displaystyle {\int }_e}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}+{\mu }^{ext}{\displaystyle {\int }_{T_e^{ext}}}{u}^{ext}\widetilde{{\sigma }_e}\end{array}$ (57)

Introduce the expressions:

$R_e\mathrm{<mo>\; :\; </mo>}=Q-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n\left(e\right)}+{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}$ (58)

$a\mathrm{<mo>\; :\; </mo>}={\displaystyle {\int }_e}{R}_e\widetilde{{\sigma }_e}={\displaystyle {\int }_e}Q\widetilde{{\sigma }_e}-{\epsilon }^{int}{\displaystyle {\int }_e}\frac{\partial u_h^{int}}{\partial n\left(e\right)}\widetilde{{\sigma }_e}+{\epsilon }^{ext}{\displaystyle {\int }_e}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}\widetilde{{\sigma }_e}$ (59)

$={\epsilon }^{int}{\displaystyle {\int }_e}\frac{\partial u^{int}}{\partial n\left(e\right)}\widetilde{{\sigma }_e}-{\epsilon }^{ext}{\displaystyle {\int }_e}\frac{\partial u^{ext}}{\partial n\left(e\right)}\widetilde{{\sigma }_e}-{\epsilon }^{int}{\displaystyle {\int }_e}\frac{\partial u_h^{int}}{\partial n\left(e\right)}\widetilde{{\sigma }_e}+{\epsilon }^{ext}{\displaystyle {\int }_e}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}\widetilde{{\sigma }_e}$ (60)

Apply now (56) and (57) to the last identity in order to obtain $a=b$ where

$\begin{array}{c}b\mathrm{<mo>\; :\; </mo>}=-{\displaystyle {\int }_{T_e^{int}}}\left(f+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)\widetilde{{\sigma }_e}-{\displaystyle {\int }_{T_e^{ext}}}\left(f+{\epsilon }^{ext}\Delta u_h^{ext}-{\mu }^{ext}{u}_h^{ext}\right)\widetilde{{\sigma }_e}\\ +{\epsilon }^{int}{\displaystyle {\int }_{T_e^{int}}}\nabla \left(u^{int}-u_h^{int}\right)\cdot \nabla \widetilde{{\sigma }_e}+{\epsilon }^{ext}{\displaystyle {\int }_{T_e^{ext}}}\nabla \left(u^{ext}-u_h^{ext}\right)\cdot \nabla \widetilde{{\sigma }_e}\\ +{\mu }^{int}{\displaystyle {\int }_{T_e^{int}}}\left(u^{int}-u_h^{int}\right)\widetilde{{\sigma }_e}+{\mu }^{ext}{\displaystyle {\int }_{T_e^{ext}}}\left(u^{ext}-u_h^{ext}\right)\widetilde{{\sigma }_e}\end{array}$ (61)

By adding ${\displaystyle {\int }_e}\left(Q_e-Q\right)\widetilde{{\sigma }_e}$ on both sides of $a=b$ , one has on the one hand

$A\mathrm{<mo>\; :\; </mo>}=a+{\displaystyle {\int }_e}\left(Q_e-Q\right)\widetilde{{\sigma }_e}={\displaystyle {\int }_e}\left(Q_e-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n\left(e\right)}+{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}\right)\widetilde{{\sigma }_e}$ (62)

$={\displaystyle {\int }_e}{\left(Q_e-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n\left(e\right)}+{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}\right)}^2{\omega }_e^{\alpha }$ (63)

$={\displaystyle {\int }_e}{\left\{\left(Q_e-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n\left(e\right)}+{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n\left(e\right)}\right){\omega }_e^{\alpha }\right\}}^2{\omega }_e^{-\alpha }$ (64)

$={\displaystyle {\int }_e}{\sigma }_e^2{\omega }_e^{-\alpha }={\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2$ (65)

On the other hand, one has

$\begin{array}{c}B=b+{\displaystyle {\int }_e}\left(Q_e-Q\right)\widetilde{{\sigma }_e}=-{\displaystyle {\int }_{T_e^{int}}}\left(f+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)\widetilde{{\sigma }_e}\\ -{\displaystyle {\int }_{T_e^{ext}}}\left(f+{\epsilon }^{ext}\Delta u_h^{ext}-{\mu }^{ext}{u}_h^{ext}\right)\widetilde{{\sigma }_e}+{\epsilon }^{int}{\displaystyle {\int }_{T_e^{int}}}\nabla \left(u^{int}-u_h^{int}\right)\cdot \nabla \widetilde{{\sigma }_e}\\ +{\epsilon }^{ext}{\displaystyle {\int }_{T_e^{ext}}}\nabla \left(u^{ext}-u_h^{ext}\right)\cdot \nabla \widetilde{{\sigma }_e}+{\mu }^{int}{\displaystyle {\int }_{T_e^{int}}}\left(u^{int}-u_h^{int}\right)\widetilde{{\sigma }_e}\\ +{\epsilon }^{ext}{\displaystyle {\int }_{T_e^{ext}}}\left(u^{ext}-u_h^{ext}\right)\widetilde{{\sigma }_e}+{\displaystyle {\int }_e}\left(Q_e-Q\right){\omega }_e^{\alpha /2}\widetilde{{\sigma }_e}{\omega }_e^{-\alpha /2}\end{array}$ (66)

Hence, we deduce

$\begin{array}{c}B\precsim {\Vert f+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\Vert }_{{\mathbb{L}}_2\left(T_e^{int}\right)}{\Vert \widetilde{{\sigma }_e}\Vert }_{{\mathbb{L}}_2\left(T_e^{int}\right)}\\ +{\Vert f+{\epsilon }^{ext}\Delta u_h^{ext}-{\mu }^{ext}{u}_h^{ext}\Vert }_{{\mathbb{L}}_2\left(T_e^{ext}\right)}{\Vert \widetilde{{\sigma }_e}\Vert }_{{\mathbb{L}}_2\left(T_e^{ext}\right)}\\ +{\left|u^{int}-u_h^{int}\right|}_{{ℍ}^1\left(T_e^{int}\right)}{\left|\widetilde{{\sigma }_e}\right|}_{{ℍ}^1\left(T_e^{int}\right)}+{\left|u^{ext}-u_h^{ext}\right|}_{{ℍ}^1\left(T_e^{ext}\right)}{\left|\widetilde{{\sigma }_e}\right|}_{{ℍ}^1\left(T_e^{ext}\right)}\\ +{\Vert u^{int}-u_h^{int}\Vert }_{{\mathbb{L}}_2\left(T_e^{int}\right)}{\Vert \widetilde{{\sigma }_e}\Vert }_{{\mathbb{L}}_2\left(T_e^{int}\right)}+{\Vert u^{ext}-u_h^{ext}\Vert }_{{\mathbb{L}}_2\left(T_e^{ext}\right)}{\Vert \widetilde{{\sigma }_e}\Vert }_{{\mathbb{L}}_2\left(T_e^{ext}\right)}\\ +{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}{\Vert \widetilde{{\sigma }_e}{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}\end{array}$ (67)

By using the definition of the extension together with the properties (50) and (51), one has

${\Vert \widetilde{{\sigma }_e}\Vert }_{{\mathbb{L}}_2\left(T_e^{int}\right)}^2\precsim h\left(T_e^{int}\right)\gamma {\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2$ (68)

${\Vert \widetilde{{\sigma }_e}\Vert }_{{\mathbb{L}}_2\left(T_e^{ext}\right)}^2\precsim h\left(T_e^{ext}\right)\gamma {\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2$ (69)

${\left|\widetilde{{\sigma }_e}\right|}_{{ℍ}^1\left(T_e^{int}\right)}^2\precsim \frac{1}{h\left(T_e^{int}\right)}\left(\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}\right){\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2$ (70)

${\left|\widetilde{{\sigma }_e}\right|}_{{ℍ}^1\left(T_e^{ext}\right)}^2\precsim \frac{1}{h\left(T_e^{ext}\right)}\left(\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}\right){\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2$ (71)

A combination of (68)-(71) with the expression of $B$ yields

$\begin{array}{c}B\precsim {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\{{\Vert f+{\epsilon }^T\Delta u_h-{\mu }^T{u}_h\Vert }_{{\mathbb{L}}_2\left(T\right)}\sqrt{\gamma h\left(T\right)}{\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}\\ +{\left|u-u_h\right|}_{{ℍ}^1\left(T\right)}\sqrt{\frac{1}{h\left(T\right)}\left(\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}\right)}{\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}^2\left(e\right)}\\ +{\Vert u-u_h\Vert }_{{\mathbb{L}}_2\left(T\right)}\sqrt{\gamma h\left(T\right)}{\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}\}\\ +{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}{\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}\end{array}$ (72)

From the equality of A and B, one obtains for the interface edge $e\in {\mathbb{E}}_h^{\Gamma }$ :

$\begin{array}{c}{\Vert {\sigma }_e{\omega }_e^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}\precsim {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\{{\Vert f+{\epsilon }^T\Delta u_h-{\mu }^T{u}_h\Vert }_{{\mathbb{L}}_2\left(T\right)}\sqrt{\gamma h\left(T\right)}\\ +{\left|u-u_h\right|}_{{ℍ}^1\left(T\right)}\sqrt{\frac{1}{h\left(T\right)}\left(\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}\right)}\\ +{\Vert u-u_h\Vert }_{{\mathbb{L}}_2\left(T\right)}\sqrt{\gamma h\left(T\right)}\}+{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}\end{array}$ (73)

Consequently, from (52) and the definition (37) of the interface estimator ${\eta }_{\alpha \mathrm{,}e}^{\Gamma }$ , one obtains for each fixed polynomial degree $p\ge 1$ :

$\begin{array}{c}\frac{2p}{h\left(e\right)}{\left({\eta }_{\alpha \mathrm{,}e}^{\Gamma }\right)}^2\precsim {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\{{\Vert f+{\epsilon }^T\Delta u_h-{\mu }^T{u}_h\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\sqrt{\gamma h\left(T\right)}\\ +{\left|u-u_h\right|}_{{ℍ}^1\left(T\right)}^2\frac{1}{h\left(T\right)}\left(\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}\right)\\ +{\Vert u-u_h\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\gamma h\left(T\right)\}+{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2\end{array}$ (74)

By inserting $f_T$ in the last right hand side and by multiplying with $h\left(e\right)/\left(2p\right)$ , one deduces

$\begin{array}{c}{\left({\eta }_{\alpha \mathrm{,}e}^{\Gamma }\right)}^2\precsim {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\{\gamma \frac{h^2\left(T\right)}{p}{\Vert f_T+{\epsilon }^T\Delta u_h-{\mu }^T{u}_h\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\\ +\left(\frac{\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}}{p}\right){\left|u-u_h\right|}_{{ℍ}^1\left(T\right)}^2+\frac{\gamma h^2\left(T\right)}{p}{\Vert f-f_T\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\\ +\frac{h^2\left(T\right)\gamma }{p}{\Vert u-u_h\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\}+\frac{h\left(e\right)}{p}{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2\end{array}$ (75)

Combine with the estimation of the local residual to obtain:

$\begin{array}{c}{\left({\eta }_{\alpha \mathrm{,}e}^{\Gamma }\right)}^2\precsim {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\{\gamma \frac{h^2\left(T\right)}{p}{\left|u-u_h\right|}_{{ℍ}^1\left(T\right)}^2\frac{p^4}{h^2\left(T\right)}+\gamma \frac{h^2\left(T\right)}{p}{\Vert f-f_T\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\\ +\left(\frac{\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}}{p}\right){\left|u-u_h\right|}_{{ℍ}^1\left(T\right)}^2+\frac{\gamma h^2\left(T\right)}{p}{\Vert f-f_T\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\\ +\frac{h^2\left(T\right)\gamma }{p}{\Vert u-u_h\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\}+\frac{h\left(e\right)}{p}{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2\end{array}$ (76)

By regrouping the terms, one obtains

$\begin{array}{c}{\left({\eta }_{\alpha \mathrm{,}e}^{\Gamma }\right)}^2\precsim {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\{\left(\gamma p^3+\frac{\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}}{p}\right){\left|u-u_h\right|}_{{ℍ}^1\left(T\right)}^2\\ +\frac{h^2\left(T\right)\gamma }{p}{\Vert u-u_h\Vert }_{{\mathbb{L}}_2\left(T\right)}^2+\frac{\gamma h^2\left(T\right)}{p}{\Vert f-f_T\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\}\\ +\frac{h\left(e\right)}{p}{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2\end{array}$ (77)

Then, for the local load oscillation ${\Vert f-f_T\Vert }_{{\mathbb{L}}_2\left(T\right)}^2$ and the local weighted flux oscillation ${\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2$ one has

$\begin{array}{c}{\eta }_{\alpha \mathrm{,}e}^{\Gamma }\precsim {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\sqrt{\gamma p^3+\frac{\gamma p^{2\left(2-\alpha \right)}+{\gamma }^{-1}}{p}+\gamma \frac{h^2\left(T\right)}{p}}{\Vert u-u_h\Vert }_{{ℍ}^1\left(T\right)}\\ +{\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\frac{\sqrt{\gamma }h\left(T\right)}{\sqrt{p}}{\Vert f-f_T\Vert }_{{\mathbb{L}}_2\left(T\right)}+\sqrt{\frac{h\left(e\right)}{p}}{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2\end{array}$ (78)

By using the assumption $h^2\le Cp$ , one concludes the last estimate in the theorem.

Theorem 2. Let $e$ be an internal edge in ${\mathbb{E}}_h^{int}$ (resp. an external edge in ${\mathbb{E}}_h^{ext}$ ). Under the assumption that $h^2\le Cp$ , we have the following bound:

$\begin{array}{c}{\eta }_{\alpha \mathrm{,}e}^{int}\le {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}{C}_{p\mathrm{,}\gamma }{\Vert u-u_h\Vert }_{{ℍ}^1\left(T\right)}+C_{\gamma }{\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\frac{h\left(T\right)}{\sqrt{p}}{\Vert f-f_T\Vert }_{{\mathbb{L}}_2\left(T\right)}\\ +\sqrt{\frac{h\left(e\right)}{p}}{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2\end{array}$ (79)

and respectively

$\begin{array}{c}{\eta }_{\alpha \mathrm{,}e}^{ext}\le {\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}{C}_{p\mathrm{,}\gamma }{\Vert u-u_h\Vert }_{{ℍ}^1\left(T\right)}+C_{\gamma }{\displaystyle \underset{T\in \mathcal{N}\left(e\right)}{\sum }}\frac{h\left(T\right)}{\sqrt{p}}{\Vert f-f_T\Vert }_{{\mathbb{L}}_2\left(T\right)}\\ +\sqrt{\frac{h\left(e\right)}{p}}{\Vert \left(Q_e-Q\right){\omega }_e^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2\end{array}$ (80)

Proof. Define

${\sigma }_e^{int}\mathrm{<mo>\; :\; </mo>}={\epsilon }^{int}〚\frac{\partial u_h^{int}}{\partial n\left(e\right)}〛{\omega }_e^{\alpha }\left(\text{resp}\mathrm{.}{\sigma }_e^{ext}\mathrm{<mo>\; :\; </mo>}={\epsilon }^{ext}〚\frac{\partial u_h^{ext}}{\partial n\left(e\right)}〛{\omega }_e^{\alpha }\right)$ (81)

and proceed as in the former proof.

Theorem 3. Let the domain oscillation be:

${\text{osc}}^{\Omega }\mathrm{<mo>\; :\; </mo>}=\sqrt{{\displaystyle \underset{T\in {\mathbb{M}}_h}{\sum }}{\Vert f-f_T\Vert }_{{\mathbb{L}}^2\left(T\right)}^2}$ (82)

and the interface oscillation be:

${\text{osc}}^{\Gamma }\mathrm{<mo>\; :\; </mo>}=\sqrt{{\displaystyle \underset{e\in {\mathbb{E}}_h^{\Gamma }}{\sum }}{\Vert Q-Q_e\Vert }_{{\mathbb{L}}^2\left(e\right)}^2}$ (83)

The next estimate holds for the weighted error estimator

${\left|u-u_h\right|}_{{ℍ}^1\left(\Omega \right)}^2\le p^{\alpha /2}\left\{{\displaystyle \underset{T\in {\mathbb{M}}_h}{\sum }}{\left({\eta }_{\alpha \mathrm{,}T}^{loc}\right)}^2\right\}+h^2{p}^2{\left({\text{osc}}^{\Omega }\right)}^2+hp{\left({\text{osc}}^{\Gamma }\right)}^2$ (84)

Proof. Due to (29), one has for $I\mathrm{<mo>\; :\; </mo>}=\mathcal{A}\left(u-u_h\mathrm{,}w\right)$ from the bilinear form (33)

$\begin{array}{c}I={\displaystyle {\int }_{\Omega }}fw+{\displaystyle {\int }_{\Gamma }}Qw-{\displaystyle \underset{T\in {\mathbb{M}}_h^{int}}{\sum }}\left\{{\displaystyle {\int }_T}\left(-{\epsilon }^{int}\Delta u_h^{int}+{\mu }^{int}{u}_h^{int}\right)w+{\epsilon }^{int}{\displaystyle {\int }_{\partial T}}\frac{\partial u_h^{ext}}{\partial n}w\right\}\\ -{\displaystyle \underset{T\in {\mathbb{M}}_h^{ext}}{\sum }}\left\{{\displaystyle {\int }_T}\left(-{\epsilon }^{ext}\Delta u_h^{ext}+{\mu }^{ext}{u}_h^{ext}\right)w+{\epsilon }^{ext}{\displaystyle {\int }_{\partial T}}\frac{\partial u_h^{ext}}{\partial n}w\right\}\end{array}$ (85)

$\begin{array}{c}I={\displaystyle {\int }_{\Gamma }}Qw+{\displaystyle \underset{T\in {\mathbb{M}}_h^{int}}{\sum }}\left\{{\displaystyle {\int }_T}fw+{\displaystyle {\int }_T}\left({\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)w-{\epsilon }^{int}{\displaystyle {\int }_{\partial T}}\frac{\partial u_h}{\partial n}w\right\}\\ +{\displaystyle \underset{T\in {\mathbb{M}}_h^{ext}}{\sum }}\left\{{\displaystyle {\int }_T}fw+{\displaystyle {\int }_T}\left({\epsilon }^{ext}\Delta v^{ext}-{\mu }^{ext}{v}^{ext}\right)w-{\epsilon }^{ext}{\displaystyle {\int }_{\partial T}}\frac{\partial u_h}{\partial n}w\right\}\end{array}$ (86)

Every side of $\partial T$ is represented as an edge of ${\mathbb{E}}_h$ which is categorized in ${\mathbb{E}}_h^{\Gamma }$ , ${\mathbb{E}}_h^0$ , ${\mathbb{E}}_h^{int}$ , ${\mathbb{E}}_h^{ext}$ . As a consequence,

$\begin{array}{c}I={\displaystyle \underset{T\in {\mathbb{M}}_h^{int}}{\sum }}\left\{{\displaystyle {\int }_T}\left(f+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)w\right\}\\ +{\displaystyle \underset{T\in {\mathbb{M}}_h^{ext}}{\sum }}\left\{{\displaystyle {\int }_T}\left(f+{\epsilon }^{ext}\Delta v^{ext}-{\mu }^{ext}{v}^{ext}\right)w\right\}+A^{int}+A^{ext}+A^0+A^{\Gamma }\end{array}$ (87)

$A^{int}\mathrm{<mo>\; :\; </mo>\; <mo>\; =\; </mo>}{\displaystyle \underset{e\in {\mathbb{E}}_h^{int}}{\sum }}{\displaystyle {\int }_e}〚{\epsilon }^{int}\frac{\partial u_h}{\partial n}〛w\mathrm{,}{A}^{ext}\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{e\in {\mathbb{E}}_h^{ext}}{\sum }}{\displaystyle {\int }_e}〚{\epsilon }^{ext}\frac{\partial u_h}{\partial n}〛w$ (88)

$A^0\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{e\in {\mathbb{E}}_h^0}{\sum }}{\displaystyle {\int }_e}\left({\epsilon }^{ext}\frac{\partial u_h}{\partial n}\right)w\mathrm{,}{A}^{\Gamma }\mathrm{<mo>\; :\; </mo>}={\displaystyle \underset{e\in {\mathbb{E}}_h^{\Gamma }}{\sum }}{\displaystyle {\int }_e}Q-\left({\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n}-{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n}\right)w$ (89)

where $A^0\equiv 0$ because ${w|}_{\partial \Omega }=0$ . Therefore,

$\begin{array}{c}I\le {\displaystyle \underset{T\in {\mathbb{M}}_h^{int}}{\sum }}\left\{{\Vert f+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\Vert }_{{\mathbb{L}}_2\left(T\right)}{\Vert w\Vert }_{{\mathbb{L}}_2\left(T\right)}\right\}\\ +{\displaystyle \underset{T\in {\mathbb{M}}_h^{ext}}{\sum }}\left\{{\Vert f+{\epsilon }^{ext}\Delta v^{ext}-{\mu }^{ext}{v}^{ext}\Vert }_{{\mathbb{L}}_2\left(T\right)}{\Vert w\Vert }_{{\mathbb{L}}_2\left(T\right)}\right\}\\ +{\displaystyle \underset{e\in {\mathbb{E}}_h^{int}}{\sum }}{\Vert 〚{\epsilon }^{int}\frac{\partial u_h}{\partial n}〛\Vert }_{{\mathbb{L}}_2\left(e\right)}{\Vert w\Vert }_{{\mathbb{L}}_2\left(e\right)}+{\displaystyle \underset{e\in {\mathbb{E}}_h^{ext}}{\sum }}{\Vert 〚{\epsilon }^{ext}\frac{\partial u_h}{\partial n}〛\Vert }_{{\mathbb{L}}_2\left(e\right)}{\Vert w\Vert }_{{\mathbb{L}}_2\left(e\right)}\\ +{\displaystyle \underset{e\in {\mathbb{E}}_h^{\Gamma }}{\sum }}{\Vert Q-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n}-{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n}\Vert }_{{\mathbb{L}}_2\left(e\right)}{\Vert w\Vert }_{{\mathbb{L}}_2\left(e\right)}\end{array}$ (90)

In particular, for $w=\left(u-u_h\right)-I_h\left(u-u_h\right)$ where $I_h$ is the Clément interpolant in [17] [19] . Since $\mathcal{A}\left(u-u_h\mathrm{,}{I}_h\left(u-u_h\right)\right)=0$ , one has $\mathcal{A}\left(u-u_h\mathrm{,}u-u_h-I_h\left(u-u_h\right)\right)=\mathcal{A}\left(u-u_h\mathrm{,}u-u_h\right)$ .

$\begin{array}{c}{\left|u-u_h\right|}_{{ℍ}^1\left(\Omega \right)}^2\le {\displaystyle \underset{T\in {\mathbb{M}}_h^{int}}{\sum }}\left\{{\Vert f+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\Vert }_{{\mathbb{L}}_2\left(T\right)}{\Vert \left(u-u_h\right)-I_h\left(u-u_h\right)\Vert }_{{\mathbb{L}}_2\left(T\right)}\right\}\\ +{\displaystyle \underset{T\in {\mathbb{M}}_h^{ext}}{\sum }}\left\{{\Vert f+{\epsilon }^{ext}\Delta v^{ext}-{\mu }^{ext}{v}^{ext}\Vert }_{{\mathbb{L}}_2\left(T\right)}{\Vert \left(u-u_h\right)-I_h\left(u-u_h\right)\Vert }_{{\mathbb{L}}_2\left(T\right)}\right\}\\ +{\displaystyle \underset{e\in {\mathbb{E}}_h^{int}}{\sum }}{\Vert 〚{\epsilon }^{int}\frac{\partial u_h}{\partial n}〛\Vert }_{{\mathbb{L}}_2\left(e\right)}{\Vert \left(u-u_h\right)-I_h\left(u-u_h\right)\Vert }_{{\mathbb{L}}_2\left(e\right)}\\ +{\displaystyle \underset{e\in {\mathbb{E}}_h^{ext}}{\sum }}{\Vert 〚{\epsilon }^{ext}\frac{\partial u_h}{\partial n}〛\Vert }_{{\mathbb{L}}_2\left(e\right)}{\Vert \left(u-u_h\right)-I_h\left(u-u_h\right)\Vert }_{{\mathbb{L}}_2\left(e\right)}\\ +{\displaystyle \underset{e\in {\mathbb{E}}_h^{\Gamma }}{\sum }}{\Vert Q-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n}+{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n}\Vert }_{{\mathbb{L}}_2\left(e\right)}{\Vert \left(u-u_h\right)-I_h\left(u-u_h\right)\Vert }_{{\mathbb{L}}_2\left(e\right)}\end{array}$ (91)

Use the interpolation property of the Clément interpolant [17] [19] to obtain

${\Vert v-I_{h}v\Vert }_{{\mathbb{L}}_2\left(T\right)}\le \frac{h\left(T\right)}{p}{\left|v\right|}_{{ℍ}^1\left(\mathcal{N}\left(T\right)\right)}\mathrm{,}{\Vert v-I_{h}v\Vert }_{{\mathbb{L}}_2\left(e\right)}\le \sqrt{\frac{h\left(e\right)}{p}}{\left|v\right|}_{{ℍ}^1\left(\mathcal{N}\left(e\right)\right)}$ (92)

Deduce therefore the next estimate

$\begin{array}{c}{\left|u-u_h\right|}_{{ℍ}^1\left(\Omega \right)}^2\le {\displaystyle \underset{T\in {\mathbb{M}}_h^{int}}{\sum }}\left\{\frac{h\left(T\right)}{p}{\Vert f+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\Vert }_{{\mathbb{L}}_2\left(T\right)}{\left|u-u_h\right|}_{{ℍ}^1\left(\mathcal{N}\left(T\right)\right)}\right\}\\ +{\displaystyle \underset{T\in {\mathbb{M}}_h^{ext}}{\sum }}\left\{\frac{h\left(T\right)}{p}{\Vert f+{\epsilon }^{ext}\Delta v^{ext}-{\mu }^{ext}{v}^{ext}\Vert }_{{\mathbb{L}}_2\left(T\right)}{\left|u-u_h\right|}_{{ℍ}^1\left(\mathcal{N}\left(T\right)\right)}\right\}\\ +{\displaystyle \underset{e\in {\mathbb{E}}_h^{int}}{\sum }}\sqrt{\frac{h\left(e\right)}{p}}{\Vert 〚{\epsilon }^{int}\frac{\partial u_h}{\partial n}〛\Vert }_{{\mathbb{L}}_2\left(e\right)}{\left|u-u_h\right|}_{{ℍ}^1\left(\mathcal{N}\left(e\right)\right)}\\ +{\displaystyle \underset{e\in {\mathbb{E}}_h^{ext}}{\sum }}\sqrt{\frac{h\left(e\right)}{p}}{\Vert 〚{\epsilon }^{ext}\frac{\partial u_h}{\partial n}〛\Vert }_{{\mathbb{L}}_2\left(e\right)}{\left|u-u_h\right|}_{{ℍ}^1\left(\mathcal{N}\left(e\right)\right)}\\ +{\displaystyle \underset{e\in {\mathbb{E}}_h^{\Gamma }}{\sum }}\sqrt{\frac{h\left(e\right)}{p}}{\Vert Q-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n}+{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n}\Vert }_{{\mathbb{L}}_2\left(e\right)}{\left|u-u_h\right|}_{{ℍ}^1\left(\mathcal{N}\left(e\right)\right)}\end{array}$ (93)

$\begin{array}{c}{\left|u-u_h\right|}_{{ℍ}^1\left(\Omega \right)}^2\le \{{\displaystyle \underset{T\in {\mathbb{M}}_h^{int}}{\sum }}\frac{h^2\left(T\right)}{p^2}{\Vert f+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\\ +{\displaystyle \underset{T\in {\mathbb{M}}_h^{ext}}{\sum }}\frac{h^2\left(T\right)}{p^2}{\Vert f+{\epsilon }^{ext}\Delta v^{ext}-{\mu }^{ext}{v}^{ext}\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\\ +{\displaystyle \underset{e\in {\mathbb{E}}_h^{int}}{\sum }}\frac{h\left(e\right)}{p}{\Vert 〚{\epsilon }^{int}\frac{\partial u_h}{\partial n}〛\Vert }_{{\mathbb{L}}_2\left(e\right)}^2+{\displaystyle \underset{e\in {\mathbb{E}}_h^{ext}}{\sum }}\frac{h\left(e\right)}{p}{\Vert 〚{\epsilon }^{ext}\frac{\partial u_h}{\partial n}〛\Vert }_{{\mathbb{L}}_2\left(e\right)}^2\\ +{{\displaystyle \underset{e\in {\mathbb{E}}_h^{\Gamma }}{\sum }}\frac{h\left(e\right)}{p}{\Vert Q-{\epsilon }^{int}\frac{\partial u_h^{int}}{\partial n}+{\epsilon }^{ext}\frac{\partial u_h^{ext}}{\partial n}\Vert }_{{\mathbb{L}}_2\left(e\right)}^2\}}^{1/2}\\ \times {\left\{{\displaystyle \underset{T\in {\mathbb{M}}_h}{\sum }}{\left|u-u_h\right|}_{{ℍ}^1\left(\mathcal{N}\left(T\right)\right)}^2\right\}}^{1/2}\end{array}$ (94)

Deduce the results for the vanishing weight $\alpha \equiv 0$ :

$\begin{array}{c}{\left|u-u_h\right|}_{{ℍ}^1\left(\Omega \right)}^2\le \{{\displaystyle \underset{T\in {\mathbb{M}}_h}{\sum }}\frac{h^2\left(T\right)}{p^2}{\Vert f-f_T\Vert }_{{\mathbb{L}}^2\left(T\right)}^2+{\displaystyle \underset{T\in {\mathbb{M}}_h^{int}}{\sum }}{\left({\eta }_{\mathrm{0,}T}^{int}\right)}^2+{\displaystyle \underset{T\in {\mathbb{M}}_h^{ext}}{\sum }}{\left({\eta }_{\mathrm{0,}T}^{ext}\right)}^2\\ +{\displaystyle \underset{e\in {\mathbb{E}}_h^{int}}{\sum }}{\left({\eta }_{\mathrm{0,}e}^{int}\right)}^2+{\displaystyle \underset{e\in {\mathbb{E}}_h^{ext}}{\sum }}{\left({\eta }_{\mathrm{0,}e}^{ext}\right)}^2+{\displaystyle \underset{e\in {\mathbb{E}}_h}{\sum }}\frac{h\left(e\right)}{p}{\Vert Q-Q_e\Vert }_{{\mathbb{L}}^2\left(e\right)}^2\\ +{{\displaystyle \underset{e\in {\mathbb{E}}_h^{int}}{\sum }}{\left({\eta }_{\mathrm{0,}e}^{\Gamma }\right)}^2\}}^{1/2}{\left|u-u_h\right|}_{{ℍ}^1\left(\Omega \right)}\end{array}$ (95)

Regroup the local terms with respect to the the local edges to obtain

${\left|u-u_h\right|}_{{ℍ}^1\left(\Omega \right)}^2\le {\displaystyle \underset{T\in {\mathbb{M}}_h}{\sum }}{\left({\eta }_{\mathrm{0,}T}^{loc}\right)}^2+{\displaystyle \underset{T\in {\mathbb{M}}_h}{\sum }}\frac{h^2\left(T\right)}{p^2}{\Vert f-f_T\Vert }_{{\mathbb{L}}^2\left(T\right)}^2+{\displaystyle \underset{e\in {\mathbb{E}}_h^{\Gamma }}{\sum }}\frac{h\left(e\right)}{p}{\Vert Q-Q_e\Vert }_{{\mathbb{L}}^2\left(e\right)}^2$ (96)

For non-vanishing weights $\alpha \ne 0$ , one applies the polynomial inverse esti- mate (47) to the polynomial

${\pi }_p\mathrm{<mo>\; :\; </mo>}=f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}$ (97)

in order to obtain

${\Vert f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\Vert }_{{\mathbb{L}}^2\left(T\right)}\le p^{2\alpha }{\Vert \left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}^2\left(T\right)}$ (98)

By using the same thing for the external domain, obtain ${\left({\eta }_{\mathrm{0,}T}^{loc}\right)}^2\le p^{2\alpha }{\left({\eta }_{\alpha \mathrm{,}T}^{loc}\right)}^2$ . Hence,

${\left|u-u_h\right|}_{{ℍ}^1\left(\Omega \right)}^2\le {\displaystyle \underset{T\in {\mathbb{M}}_h}{\sum }}{p}^{\alpha /2}{\left({\eta }_{\alpha \mathrm{,}T}^{loc}\right)}^2+{\displaystyle \underset{T\in {\mathbb{M}}_h}{\sum }}\frac{h^2\left(T\right)}{p^2}{\Vert f-f_T\Vert }_{{\mathbb{L}}^2\left(T\right)}^2+{\displaystyle \underset{e\in {\mathbb{E}}_h^{\Gamma }}{\sum }}\frac{h\left(e\right)}{p}{\Vert Q-Q_e\Vert }_{{\mathbb{L}}^2\left(e\right)}^2$ (99)

Theorem 4. For an element $T\in {\mathbb{M}}_h^{int}$ , respectively $T\in {\mathbb{M}}_h^{ext}$ , introduce:

${\sigma }_T^{int}\mathrm{<mo>\; :\; </mo>}=\left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right){\omega }_T^{\alpha }\mathrm{,}{\sigma }_T^{ext}\mathrm{<mo>\; :\; </mo>}=\left(f_T+{\epsilon }^{ext}\Delta u_h^{ext}-{\mu }^{ext}{u}_h^{ext}\right){\omega }_T^{\alpha }$ (100)

One has for a fixed polynomial degree $p\ge 1$ the following estimates

${\Vert {\sigma }_T^{int}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}\precsim \frac{p^{\left(2-\alpha \right)}}{h\left(T\right)}{\Vert u^{int}-u_h^{int}\Vert }_{{ℍ}_1\left(T\right)}+{\Vert \left(f_T-f\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}$ (101)

${\Vert {\sigma }_T^{ext}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}\precsim \frac{p^{\left(2-\alpha \right)}}{h\left(T\right)}{\Vert u^{ext}-u_h^{ext}\Vert }_{{ℍ}_1\left(T\right)}+{\Vert \left(f_T-f\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}$ (102)

Thus, the local expressions ${\eta }_{\alpha \mathrm{,}T}^{int}$ and ${\eta }_{\alpha \mathrm{,}T}^{ext}$ verify:

${\eta }_{\alpha \mathrm{,}T}^{int}\le C_{p\mathrm{,}\alpha }{\Vert u^{int}-u_h^{int}\Vert }_{{ℍ}_1\left(T\right)}+\frac{h\left(T\right)}{p}{\Vert \left(f_T-f\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}$ (103)

${\eta }_{\alpha \mathrm{,}T}^{ext}\le C_{p\mathrm{,}\alpha }{\Vert u^{ext}-u_h^{ext}\Vert }_{{ℍ}_1\left(T\right)}+\frac{h\left(T\right)}{p}{\Vert \left(f_T-f\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}$ (104)

Proof. The following equalities hold:

${\Vert {\sigma }_T^{int}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}={\displaystyle {\int }_T}{\left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)}^2{\omega }_T^{\alpha }$ (105)

$={\displaystyle {\int }_T}\left(f+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right){\sigma }_T^{int}+{\displaystyle {\int }_T}\left(f_T-f\right){\sigma }_T^{int}$ (106)

$\begin{array}{l}={\epsilon }^{int}{\displaystyle {\int }_T}\nabla u^{int}\cdot \nabla {\sigma }_T^{int}+{\mu }^{int}{\displaystyle {\int }_T}{u}^{int}{\sigma }_T^{int}-{\epsilon }^{int}{\displaystyle {\int }_T}\nabla u_h^{int}\cdot \nabla {\sigma }_T^{int}\\ -{\mu }^{int}{\displaystyle {\int }_T}{u}_h^{int}{\sigma }_T^{int}+{\displaystyle {\int }_T}\left(f_T-f\right){\sigma }_T^{int}\end{array}$ (107)

$\begin{array}{l}={\epsilon }^{int}{\displaystyle {\int }_T}\nabla \left(u^{int}-u_h^{int}\right)\cdot \nabla {\sigma }_T^{int}+{\mu }^{int}{\displaystyle {\int }_T}{\left(u-u_h\right)}^{int}{\sigma }_T^{int}\\ +{\displaystyle {\int }_T}\left(f_T-f\right){\omega }_T^{\alpha /2}{\sigma }_T^{int}{\omega }_T^{-\alpha /2}\end{array}$ (108)

Consequently, one obtains

${\Vert {\sigma }_T^{int}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}\precsim {\Vert u^{int}-u_h^{int}\Vert }_{{ℍ}^1\left(T\right)}{\left|{\sigma }_T\right|}_{{ℍ}^1\left(T\right)}+{\Vert \left(f_T-f\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}{\Vert {\sigma }_T^{int}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}$ (109)

Concerning the estimation of ${\left|{\sigma }_T^{int}\right|}_{{ℍ}^1\left(T\right)}$ , one has

${\left|{\sigma }_T^{int}\right|}_{{ℍ}^1\left(T\right)}={\displaystyle {\int }_T}{\left|\nabla \left\{\left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right){\omega }_T^{\alpha }\right\}\right|}^2$ (110)

$\begin{array}{l}\le {\displaystyle {\int }_T}{\omega }_T^{2\alpha }{\left|\nabla \left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)\right|}^2\\ +{\displaystyle {\int }_T}{\left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)}^2{\left|\nabla {\omega }_T^{2\alpha }\right|}^2\end{array}$ (111)

For the first term, apply the polynomial inverse estimate (48) to the polynomial

${\pi }_p\mathrm{<mo>\; :\; </mo>}=f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}$ (112)

and the affine transform $F_T\mathrm{<mo>\; :\; </mo>}\hat{T}\to T$ in order to obtain

${\displaystyle {\int }_T}{\omega }_T^{2\alpha }{\left|\nabla \left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)\right|}^2\precsim \frac{p^{2\left(2-\alpha \right)}}{h^2\left(T\right)}{\displaystyle {\int }_T}{\omega }_T^{\alpha }{\left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)}^2$ (113)

As for the second term, split ${\left|\nabla {\omega }_T^{\alpha }\right|}^2$ into ${\left({\partial }_x{\omega }_T^{\alpha }\right)}^2$ and ${\left({\partial }_y{\omega }_T^{\alpha }\right)}^2$ , use the boundedness of ${\omega }_T^{\alpha }$ and apply the inverse estimate (47) to the polynomial (112) to obtain

${\displaystyle {\int }_T}{\left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)}^2{\left|\nabla {\omega }_T^{2\alpha }\right|}^2\precsim \frac{p^{2\left(2-\alpha \right)}}{h^2\left(T\right)}{\displaystyle {\int }_T}{\left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)}^2{\omega }_T^{\alpha }$ (114)

A combination of (111), (113) and (114) yields

${\left|{\sigma }_T^{int}\right|}_{{ℍ}^1\left(T\right)}\precsim \frac{p^{2\left(2-\alpha \right)}}{h^2\left(T\right)}{\displaystyle {\int }_T}{\left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right)}^2{\omega }_T^{\alpha }$ (115)

$\begin{array}{l}=\frac{p^{2\left(2-\alpha \right)}}{h^2\left(T\right)}{\displaystyle {\int }_T}{\left\{\left(f_T+{\epsilon }^{int}\Delta u_h^{int}-{\mu }^{int}{u}_h^{int}\right){\omega }_T^{\alpha }\right\}}^2{\omega }_T^{-\alpha }\\ =\frac{p^{2\left(2-\alpha \right)}}{h^2\left(T\right)}{\Vert {\sigma }_T^{int}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\end{array}$ (116)

As a consequence, one has

$\begin{array}{c}{\Vert {\sigma }_T^{int}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\precsim \frac{p^{\left(2-\alpha \right)}}{h\left(T\right)}{\Vert u^{int}-u_h^{int}\Vert }_{{ℍ}^1\left(T\right)}{\Vert {\sigma }_T^{int}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}\\ +{\Vert \left(f_T-f\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}{\Vert {\sigma }_T^{int}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}\end{array}$ (117)

from which (101) follows. Thus, one deduces

${\Vert {\sigma }_T^{int}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}^2\precsim \frac{p^{\left(2-\alpha \right)}}{h\left(T\right)}{\Vert u^{int}-u_h^{int}\Vert }_{{ℍ}_1\left(T\right)}+{\Vert \left(f_T-f\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}$ (118)

By using the definition (35) of the interior estimator, one has

${\left({\eta }_{\alpha \mathrm{,}T}^{int}\right)}^2=\frac{h^2\left(T\right)}{p^2}{\Vert {\sigma }_T^{int}{\omega }_T^{-\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}^2$ (119)

$\precsim \frac{p^{2\left(2-\alpha \right)}}{p^2}{\Vert u^{int}-u_h^{int}\Vert }_{{ℍ}_1\left(T\right)}+\frac{h^2\left(T\right)}{p^2}{\Vert \left(f_T-f\right){\omega }_T^{\alpha /2}\Vert }_{{\mathbb{L}}_2\left(T\right)}$ (120)

The other estimates are obtained in a similar manner.

4. Practical Results

The computer implementation is performed by a combination of C functions and C++ classes. Some LAPACK and BLAS routines are used sometimes to perform various linear algebraic operations.

4.1. Exact Precision

In this subsection, we concentrate fully on the exact precision for the purpose of obtaining insight and confidence about the accuracy of the results of the computer implementation. That is, we do not consider yet any description of the error estimator ${\eta }_{\alpha }$ . We examine several parameters comprising the polynomial degree and the problem coefficients $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]$ . The ε-ratio and the μ-ratio are the positive constants $min\left\{{\epsilon }^{int}/{\epsilon }^{ext}\mathrm{,}{\epsilon }^{ext}/{\epsilon }^{int}\right\}$ and $min\left\{{\mu }^{int}/{\mu }^{ext}\mathrm{,}{\mu }^{ext}/{\mu }^{int}\right\}$ respectively. Since the FEM-level is used extensively, we introduce it very rapidly. For a given mesh ${\mathbb{M}}_h$ on the current level $\mathcal{l}$ , another mesh on the next level $\left(\mathcal{l}+1\right)$ is constructed by subdividing every edge in the middle of which one inserts a new node. Therefore, that corresponds to a global uniform refinement where every element is locally subdivided. Only the mesh on the coarsest level (level $\mathcal{l}=1$ in our entire study) is provided. Thus, one level incrementation amounts to reducing the mesh size from h to h/2. For the numerical compu- tations, we consider the unit square ${\left[\mathrm{0,1}\right]}^2$ as the entire domain $\Omega$ while the internal domain ${\Omega }^{int}$ is ${\left[1/3\mathrm{,}2/3\right]}^2$ . The used mesh ${\mathbb{M}}_h$ is a tensor product uniform triangular mesh. The exact solution is chosen to be $\left(1/{\epsilon }^{int}\right)sin\left(3\text{π}x\right)sin\left(3\text{π}y\right)$ for the internal domain while it is $\left(1/{\epsilon }^{ext}\right)sin\left(3\text{π}x\right)sin\left(3\text{π}y\right)$ for the external one. It is globally continuous and it admits highly discontinuous derivatives at the interface $\Gamma =\partial {\Omega }^{int}$ depending on the values of ${\epsilon }^{int}$ and ${\epsilon }^{ext}$ . The right hand side function is obtained by applying the transmission operators (2) and (3) to the exact solution. The ${\mathbb{L}}_{\infty }$ -norm is practically obtained by considering a very dense point sample $\left\{x_i^T\right\}$ inside every element $T\in {\mathbb{M}}_h$ . The maximum values of $\left|u\left(x_i^T\right)-u_h\left(x_i^T\right)\right|$ over all samples and over all elements will then be the practical ${\mathbb{L}}_{\infty }$ -norm.

The first configuration for the investigation of the exact precision corresponds to $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{1000,1,0.1,100}\right]$ which associates to both the ε-ratio and the μ-ratio the value of $\mathrm{1\; <mo>\; :\; </mo>1000}$ . We use these parameters because they highlight a situation where the internal and the external coefficients are very dissimilar. In particular, the ratio of the normal derivatives at the interface is proportional to the ε-ratio. Parameters tending to unity are very easy because that turns out to be the same as treating a problem without an interface in the whole domain $\Omega$ . The results of the computation are collected in Table 1 where we consider three fixed polynomial degrees $p=1,2,3$ . For each degree, the FEM-levels are allowed to vary from one to five. The corresponding ${\mathbb{L}}_2$ - norms, ${ℍ}^1$ -seminorms and ${\mathbb{L}}_{\infty }$ -norms are depicted there as well as the contraction which is the ratio between two consecutive errors as the FEM-level is incremented. The ${ℍ}^1$ -seminorm errors are in general two digit worse than the ${\mathbb{L}}_2$ -norm. That is observed for different configurations related to the problem coefficients $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]$ and different simulation parameters. The ${\mathbb{L}}_{\infty }$ - error is just a little worse than the ${\mathbb{L}}_2$ -error but not as imprecise as the error in the ${ℍ}^1$ -seminorm. Since the ${ℍ}^1$ -norm is quadratically the sum of the ${\mathbb{L}}_2$ - norm and the ${ℍ}^1$ -seminorm as ${\left|\cdot \right|}_{{ℍ}^1\left(\Omega \right)}^2={\Vert \cdot \Vert }_{{\mathbb{L}}_2\left(\Omega \right)}^2+{\left|\cdot \right|}_{{ℍ}^1\left(\Omega \right)}^2$ , the ${ℍ}^1$ -norm and the ${ℍ}^1$ -seminorm are practically of the same order because the ${\mathbb{L}}_2$ -norm is dominated by the ${ℍ}^1$ -seminorm. In the next description, we will be mainly interested in the ${ℍ}^1$ -norm or practically the ${ℍ}^1$ -seminorm. We observe practically constant contraction numbers depending on the used polynomial degree. The same test has been conducted for another configuration of $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{0.1,100,6.5,0.5}\right]$ whose ε-ratio and μ-ratio are respec- tively $\mathrm{1\; <mo>\; :\; </mo>1000}$ and $\mathrm{1\; <mo>\; :\; </mo>13}$ . The outcome of the errors is collected in Table 2. In

Table 1. Exact precision and contraction ratio for $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{1000,1,0.1,100}\right]$ .

Table 2. Exact precision and contraction ratio for $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{0.1,100,6.5,0.5}\right]$ .

contrast to the previous case, the value of ${\epsilon }^{int}$ is currently smaller that ${\epsilon }^{ext}$ , although the ε-ratio remains the same. Nonetheless, the general characteristic of the errors remains unchanged which demonstrates the robustness of the method when put in practice.

Our next simulation about a-priori error estimates focuses on the influence of the polynomial degree p. For that, we consider four configurations corresponding to $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{0.5,2,4,0.1}\right]$ , $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{1.5,500,0.5,7}\right]$ , $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{2,0.5,1,4}\right]$ , $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{500,1.5,7,0.5}\right]$ . They have respectively the following e-ratios $\mathrm{1\; <mo>\; :\; </mo>4}$ , $\mathrm{1\; <mo>\; :\; </mo>333.3}$ , $\mathrm{1\; <mo>\; :\; </mo>4}$ , $\mathrm{1\; <mo>\; :\; </mo>333.3}$ while the μ-ratios are $\mathrm{1\; <mo>\; :\; </mo>40}$ , $\mathrm{1\; <mo>\; :\; </mo>14}$ , $\mathrm{1\; <mo>\; :\; </mo>4}$ , $\mathrm{1\; <mo>\; :\; </mo>14}$ . The corresponding ${\mathbb{L}}_2$ -errors are depicted in Figure 2(b) while the errors using the ${ℍ}^1$ -seminorm are displayed in Figure 2(a). This situation confirms again the fact that the ${\mathbb{L}}_2$ -errors are about two digit worse than the ${ℍ}^1$ -errors. In those figures, the errors are computed in function of increasing polynomial degrees which range from one to six. We observe satisfactory error decrease for all considered four configurations. Both the ${\mathbb{L}}_2$ -error and the ${ℍ}_1$ -error decrease in exponential rate with respect to the polynomial degrees. In fact, the ${\mathbb{L}}_2$ -errors drop from an order of 0.1 to an order of $1.0e-07$ in the range of $p=1,\cdots ,6$ while the ${ℍ}^1$ -errors drop from an order of five to an order of $1.0e-5$ as the polynomial degree grows.

(a)(b)

Figure 2. A-priori error in function of the increasing polynomial degree for various values of $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]$ : (a) ${\mathbb{L}}_2$ -error, (b) ${ℍ}_1$ -error.

4.2. Practical A-Posteriori Error Estimation

Now that we have gained insight about the a-priori accuracy, we want to turn our attention to the a-posteriori estimator ${\eta }_{\alpha }$ . The principal role of an a-posteriori error estimator is twofold. For one, it serves as gaining some idea of whether to continue or to abort a simulation. The computation is aborted when the desired precision is provided by the estimator. Since the exact solution u is not known for real applications, an estimator is needed. A further purpose of the error estimator is to identify the regions within the domain $\Omega$ where the precision is unsatisfactory. Mesh refinements are therefore applied at those regions to improve the local precision. The proposed estimator can be used for both purposes. But we have currently only the first utility in our computer implementation. We use different values of $\alpha$ which comprise $\alpha =0$ , $\alpha =0.25$ , $\alpha =0.5$ , $\alpha =0.75$ , $\alpha =1$ . The results of the computer simulation is collected in Table 3 where we examine four polynomial degrees $p=1,2,3,4$ . For every fixed polynomial degree p, the value of the FEM-level is allowed to

Table 3. Exact and estimated accuracies for $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{10,1,5,2}\right]$ .

range from one to five. We consider a configuration where the problem coefficients are $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{10,1,5,2}\right]$ that corresponds to an e-ratio and μ-ratio of $\mathrm{1\; <mo>\; :\; </mo>10}$ and $\mathrm{1\; <mo>\; :\; </mo>2.5}$ . According to the proposed theory, it is sufficient to conduct a comparison between the ${ℍ}^1$ -error and the error estimator ${\eta }_{\alpha }$ . For all values of $\alpha$ , Table 3 highlights that the estimator ${\eta }_{\alpha }$ provides an efficient estimation of the exact precision as predicted theoretically. In addition, the decrease of the exact precision agrees well with the decrease of the error estimator ${\eta }_{\alpha }$ as the FEM-levels grow. That can be very well observed for various values of the fixed polynomial degree p. In fact, the estimations by ${\eta }_{\alpha }$ are somewhat influenced by the values of the chosen $\alpha$ . Nonetheless, the values of ${\eta }_{\alpha }$ are comparatively of the same order up to some constant scaling factors. The same tests but for other configurations yield the results in Table 4 where we have $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{0.1,100,6.5,0.5}\right]$ . Comparable charac- teristics as in the preceding investigation are observed. That reveals again that the estimator is robust with respect to the problem configurations.

Table 4. Exact and estimated accuracies for $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{0.1,100,6.5,0.5}\right]$ .

At this point, we do no know exactly which value of the parameter $\alpha$ to be chosen. There is no mathematical study to adjust it, let alone to find its optimal value. It is conceivable to choose the value of $\alpha$ adaptively but we do not have that utility at the time being. As shown formerly, all values of $\alpha$ in $\left[\mathrm{0,1}\right]$ provide an efficient a-posteriori error estimator. Thus, the purpose of the next study is to examine the comparison of the exact accuracy in term of the ${ℍ}^1$ - norm on the the one hand with the average of the estimators ${\eta }_{\alpha }$ on the other

hand. In our case, the average is expressed as ${\eta }_{\text{avr}}=\left({\displaystyle {\sum }_{i=1}^N}{\eta }_{{\alpha }_i}\right)/N$ where we set

$N=5$ while the values of $\alpha$ are $\left\{\mathrm{0,0.25,0.5,0.75,1}\right\}$ . We examine in Figure 3(a) and Figure 3(b) the agreements of the exact precision ${\Vert u-u_h\Vert }_{{ℍ}^1\left(\Omega \right)}$ and estimated precision ${\eta }_{\text{avr}}$ when the FEM-level increases. In fact, we consider five levels for each test. Each curve corresponds to a fixed polynomial degree $p=1,2,3,4$ . The problem parameters for the coefficients in the transmission PDE are $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{0.1,100,6.5,0.5}\right]$ for the first test which is displayed in Figure 3(a). Those coefficients correspond to $\mathrm{1\; <mo>\; :\; </mo>1000}$ and $\mathrm{1\; <mo>\; :\; </mo>13}$

(a)(b)

Figure 3. Comparison of the exact accuracy and the average ${\eta }_{\text{avr}}$ in function of the FEM-levels: (a) $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{0.1,100,6.5,0.5}\right]$ ; (b) $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{10,1,5,2}\right]$ .

as ε-ratio and μ-ratio. As for the second test, we consider the problem coefficients $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]=\left[\mathrm{10,1,5,2}\right]$ which yield the ε-ratio and the μ-ratio of $\mathrm{1\; <mo>\; :\; </mo>10}$ and $\mathrm{1\; <mo>\; :\; </mo>2.5}$ respectively. An observation is that for both configurations, the curves for the exact accuracies and those for the estimated ones almost agree in shape and slope. That means that the exact accuracy is comparable to the estimated precision up to a constant factor. This highlights that using the average estimator ${\eta }_{\text{avr}}$ makes sense to evaluate the desired precision. The constant factors are problem dependent in our case as they are affected by the problem parameters $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]$ . As it can be observed in Figure 3(a), the estimated accuracies are in general somewhat smaller than the exact ones. In contrast, for the second test in Figure 3(b) where only the coefficients $\left[{\epsilon }^{int}\mathrm{,}{\epsilon }^{ext}\mathrm{,}{\mu }^{int}\mathrm{,}{\mu }^{ext}\right]$ differ from the first test, the estimated precisions are in general somewhat above the exact precisions.

5. Conclusion

We considered an interface problem where the internal and external PDEs admit different coefficients. The solution to the PDE is globally continuous but it may admit highly discontinuous derivatives across the interface separating the internal and the external domains. We proposed an a-posteriori error estimator which has been theoretically demonstrated to be equivalent to the exact precision. The performance of the estimator has also been investigated practically. For every fixed polynomial degree, it provides satisfactory precisions which are comparable to the order of the exact accuracies on different FEM-levels. Several tests have been performed where one varies the problem configurations and the simulation parameters.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Cancès, E., Mennucci, B. and Tomasi, J. (1998) A New Integral Equation Formalism for the Polarizable Continuum Model, Theoretical Background and Applications to Isotropic and Anisotropic Dielectrics. The Journal of Chemical Physics, 107, 3032-3041. https://doi.org/10.1063/1.474659
[2]	Teso, A., Filho, E. and Neto, A. (1997) Solution of the Poisson-Boltzmann Equation for a System with Four Ionic Species. Journal of Mathematical Biology, 35, 814-824. https://doi.org/10.1007/s002850050078
[3]	Tomasi, J., Mennucci, B. and Cammi, R. (2005) Quantum Mechanical Continuum Solvation Models. Chemical Reviews, 105, 2999-3094. https://doi.org/10.1021/cr9904009
[4]	Xie, D. and Zhou, S. (2007) A New Minimization Protocol for Solving Nonlinear Poisson-Boltzmann Mortar Finite Element Equation. BIT Numerical Mathematics, 47, 853-871. https://doi.org/10.1007/s10543-007-0145-9
[5]	Verfürth, R. (1994) A-Posteriori Error Estimation and Adaptive Mesh Refinement Techniques. Journal of Computational and Applied Mathematics, 50, 67-83. https://doi.org/10.1016/0377-0427(94)90290-9
[6]	Verfürth, R. (1989) A-Posteriori Error Estimators for the Stokes Equations. Numerische Mathematik, 55, 309-325. https://doi.org/10.1007/BF01390056
[7]	Verfürth, R. (1991) A-Posteriori Error Estimators for the Stokes Equations II, Non-Conforming Discretizations. Numerische Mathematik, 60, 235-249. https://doi.org/10.1007/BF01385723
[8]	Ainsworth, M. and Oden, T. (1997) A-Posteriori Error Estimators for the Stokes and Oseen Equations. SIAM Journal on Numerical Analysis, 34, 228-245. https://doi.org/10.1137/S0036142994264092
[9]	Jou, J. and Liu, J. (2000) An A-Posteriori Finite Element Analysis for the Stokes Equations. Journal of Computational and Applied Mathematics, 114, 333-343. https://doi.org/10.1016/S0377-0427(99)00271-X
[10]	Kunert, G. and Verfürth, R. (2000) Edge Residuals Dominate A-Posteriori Error Estimates for Linear Finite Element Methods on Anisotropic Triangular and Tetrahedral Meshes. Numerische Mathematik, 86, 283-303. https://doi.org/10.1007/PL00005407
[11]	Dobrowolski, M., Graf, S. and Pflaum, C. (1999) A-Posteriori Error Estimators in the Finite Element Method on Anisotropic Meshes. ETNA, 8, 36-45.
[12]	Randrianarivony, M. (2001) Stability of Mixed Finite Element Methods on Anisotropic Meshes. Master’s Thesis, Faculty of mathematics, Technische Universitat Chemnitz.
[13]	Kunert, G. (2001) Robust a Posteriori Error Estimation for a Singularly Perturbed Reaction-Diffusion Equation on Anisotropic Tetrahedral Meshes. Advances in Computational Mathematics, 15, 237-259.
[14]	Kunert, G. (2001) A-Posteriori L2 Error Estimation on Anisotropic Tetrahedral Finite Element Meshes. IMA Journal of Numerical Analysis, 21, 503-523. https://doi.org/10.1093/imanum/21.2.503
[15]	Creusé, E., Kunert, G. and Nicaise, S. (2004) A-Posteriori Error Estimation for the Stokes Problem, Anisotropic and Isotropic Discretizations. Mathematical Models and Methods in Applied Sciences, 14, 1297-1341. https://doi.org/10.1142/S0218202504003635
[16]	Siebert, K. (1996) An a Posteriori Error Estimator for Anisotropic Refinement. Numerische Mathematik, 73, 373-398. https://doi.org/10.1007/s002110050197
[17]	Melenk, J. and Wohlmuth, B. (2002) On Residual-Based a Posteriori Error Estimation in hp-FEM. Advances in Computational Mathematics, 15, 311-331.
[18]	Bernardi, C. (1996) Indicateurs d’Erreur en h-N Version des éléments Spectraux. M2AN, 30, 1-38. https://doi.org/10.1051/m2an/1996300100011
[19]	Melenk, J. (2005) hp-Interpolation of Non-Smooth Functions. SIAM Journal on Numerical Analysis, 43, 127-155. https://doi.org/10.1137/S0036142903432930
[20]	Randrianarivony, M. (2004) Anisotropic Finite Elements for the Stokes Problem, A-Posteriori Error Estimator and Adaptive Mesh. Journal of Computational and Applied Mathematics, 169, 255-275. https://doi.org/10.1016/j.cam.2003.12.025
[21]	Costabel, M. and Stephan, E. (1985) A Direct Boundary Integral Equation Method for Transmission Problems. Journal of Mathematical Analysis and Applications, 106, 367-413. https://doi.org/10.1016/0022-247X(85)90118-0
[22]	Harbrecht, H. and Randrianarivony, M. (2009) Wavelet BEM on Molecular Surfaces, Parametrization and Implementation. Computing, 86, 1-22. https://doi.org/10.1007/s00607-009-0050-y
[23]	Randrianarivony, M. (2016) Multiwavelet Boundary Element Method for Cavities Admitting Many NURBS Patches. Modeling and Numerical Simulation, 6, 69-93. https://doi.org/10.4236/mnsms.2016.64007
[24]	Harbrecht, H. and Randrianarivony, M. (2011) Wavelet BEM on Molecular Surfaces, Solvent Exclusive Surfaces. Computing, 92, 335-364. https://doi.org/10.1007/s00607-011-0147-y
[25]	Weijo, V., Randrianarivony, M., Harbrecht, H. and Frediani, L. (2010) Wavelet Formulation of the Polarizable Continuum Model. Journal of Computational Chemistry, 31, 1469-1477.
[26]	Ying, W. (2017) A Kernel-Free Boundary Integral Method for the Nonlinear Poisson-Boltzmann Equation. Journal of Computational Physics.
[27]	Holst, M., McCammon, J., Yu, Z., Zhou, Y. and Zhu, Y. (2012) Adaptive Finite Element Modeling Techniques for the Poisson-Boltzmann Equation. Communications in Computational Physics, 11, 179-214. https://doi.org/10.4208/cicp.081009.130611a
[28]	Lu, B., Zhou, Y., Holst, M. and McCammon, J. (2008) Recent Progress in Numerical Methods for the Poisson-Boltzmann Equation in Biophysical Applications. Computer Physics Communications, 3, 973-1009.
[29]	Baker, N., Sept, D., Holst, M. and McCammon, J. (2001) The Adaptive Multilevel Finite Element Solution of the Poisson-Boltzmann Equation on Massively Parallel Computers. IBM Journal of Research and Development, 45, 427-438. https://doi.org/10.1147/rd.453.0427
[30]	Chen, L., Holst, M. and Xu, J. (2007) The Finite Element Approximation of the Nonlinear Poisson-Boltzmann Equation. SIAM Journal on Numerical Analysis, 45, 2298-2320. https://doi.org/10.1137/060675514
[31]	Rabenstein, B. (2000) Monte Carlo Methods for Simulation of Protein Folding and Titration. PhD Thesis, Freie Universtat, Berlin.
[32]	Randrianarivony, M. (2006) Geometric Processing of CAD Data and Meshes as Input of Integral Equation Solvers. PhD Thesis, Technical University of Chemnitz.
[33]	Randrianarivony, M. (2009) On Global Continuity of Coons Mappings in Patching CAD Surfaces. Computer-Aided Design, 41, 782-791. https://doi.org/10.1016/j.cad.2009.04.012
[34]	Randrianarivony, M. and Brunnett, G. (2008) Preparation of CAD and Molecular Surfaces for Meshfree Solvers. Lecture Notes in Computational Science and Engineering, 65, 231-245. https://doi.org/10.1007/978-3-540-79994-8_14
[35]	Randrianarivony, M. (2008) Harmonic Variation of Edge Size in Meshing CAD Geometries from IGES Format. Lecture Notes in Computer Science, 5102, 56-65. https://doi.org/10.1007/978-3-540-69387-1_7
[36]	Randrianarivony, M. (2016) Domain Decomposition for Wavelet Single Layer on Geometries with Patches. Applied Mathematics, 7, 1798-1823. https://doi.org/10.4236/am.2016.715151
[37]	Randrianarivony, M. (2017) Analytical Polarizable Continuum Model for Wavelets on NURBS Patches. Applied Mathematics, 8, 1045-1073. https://doi.org/10.4236/am.2017.88081
[38]	Diedrich, C., Dijkstra, D., Hamaekers, J., Henninger, B. and Randrianarivony, M. (2015) A Finite Element Study on the Effect of Curvature on the Reinforcement of Matrices by Randomly Distributed and Curved Nanotubes. Journal of Computational and Theoretical Nanoscience, 12, 2108-2116. https://doi.org/10.1166/jctn.2015.3995
[39]	Wohlmuth, B. (1999) A Residual Based Error Estimator for Mortar Finite Element Discretization. Numerische Mathematik, 84, 143-171. https://doi.org/10.1007/s002110050467
[40]	Maday, Y. (1990) Analysis of Spectral Projectors in One-Dimensional Domains. Mathematics of Computation, 55, 537-562. https://doi.org/10.1090/S0025-5718-1990-1035939-1
[41]	Bernardi, C. and Maday, Y. (1992) Polynomial Interpolation Results in Sobolev Spaces. Journal of Computational and Applied Mathematics, 43, 53-80. https://doi.org/10.1016/0377-0427(92)90259-Z
[42]	Bernardi, C., Dauge, M. and Maday, Y. (2003) Polynomials in the Sobolev World. Preprint of the Laboratoire Jacques-Louis Lions, No. R03038.