Adaptive reconstruction for electrical impedance tomography with a piecewise constant conductivity

Electrical impedance tomography (EIT) aims at recovering the electrical conductivity distribution of an object from voltage measurements on the boundary. It has attracted much interest in medical imaging, geophysical prospecting, nondestructive evaluation and pneumatic oil pipeline conveying, etc. A large number of reconstruction algorithms have been proposed; see, e.g. [3, 15, 18, 21, 24, 26, 28, 30, 31, 33–37, 39, 40, 51, 54, 57] for a rather incomplete list. One prominent idea underpinning many imaging algorithms is regularization, especially Tikhonov regularization [29]. In practice, they are customarily implemented using the continuous piecewise linear finite element method (FEM) on quasi-uniform meshes, due to its flexibility in handling spatially variable coefficients and general domain geometry. The convergence analysis of finite element approximations was carried out in [22, 27, 47].

In several practical applications, the physical process is accurately described by the complete electrode model (CEM) [14, 50]. It employs nonstandard boundary conditions to capture characteristics of the experiment. In particular, around the electrode edges, the boundary condition changes from the Neumann to Robin type, which, according to classical elliptic regularity theory [23], induces weak solution singularity around the electrode edges; see, e.g. [45] for an early study. Further, the low regularity of the sought-after conductivity distribution, especially that enforced by a nonsmooth penalty, e.g. total variation, can also induce weak interior singularities of the solution. Thus, a (quasi-)uniform triangulation of the domain can be inefficient for resolving these singularities, and the discretization errors around electrode edges and internal interfaces can potentially compromise the reconstruction accuracy. These observations motivate the use of an adaptive strategy to achieve the desired accuracy in order to enhance the overall computational efficiency.

For direct problems, the mathematical theory of adaptive FEM (AFEM), including a posteriori error estimation, convergence and computational complexity, has advanced greatly [1, 12, 44, 53]. A common AFEM consists of the following successive loops:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{afem_loop} {{\rm SOLVE}}\rightarrow{{\rm ESTIMATE}}\rightarrow{{\rm MARK}}\rightarrow{{\rm REFINE}}. \nonumber \end{align} \tag{ 1.1 }$$

The module ESTIMATE employs the given problem data and computed solutions to provide computable quantities on the local errors, and distinguishes different adaptive algorithms.

In this work, we develop an adaptive EIT reconstruction algorithm with a piecewise constant conductivity. In practice, the piecewise constant nature is commonly enforced by a total variation (TV) penalty. However, it is challenging for AFEM treatment (see e.g. [5] for image denoising). Thus, we take an indirect approach based on a Modica–Mortola-type functional:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\eps}{\varepsilon} \newcommand{\Om}{\Omega} \newcommand{\dx}{\,{\rm d}x} \displaystyle \mathcal{F}_\varepsilon(\sigma) = \eps\int_\Om |\nabla \sigma|^2 \dx + \frac{1}{\eps}\int_{\Om} W(\sigma) \dx, \nonumber \end{align*}$$

where the constant $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \eps>0$ is small and $W(s):\mathbb{R}\rightarrow\mathbb{R}$ is the double-well potential, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:double-well} W(s) = (s-c_0)^2(s-c_1)^2, \nonumber \end{align} \tag{ 1.2 }$$

with $c_0,c_1>0$ being two known values that the conductivity $\sigma$ can take. The functional $\mathcal{F}_{\varepsilon}$ $\Gamma$ converges to the TV semi-norm [2, 42, 43]. The corresponding regularized least-squares formulation reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\J}{\mathcal{J}} \displaystyle \label{eqn:tikh-MM-0} \inf_{\sigma\in \widetilde{\mathcal{A}}} \left\{\J_{\varepsilon}(\sigma) = \frac{1}{2} \|U(\sigma)-U^\delta\|^2 + \frac{\widetilde{\alpha}}{2}\mathcal{F}_\varepsilon(\sigma)\right\}, \nonumber \end{align} \tag{ 1.3 }$$

where $\tilde \alpha>0$ is a regularization parameter; see section 2 for further details. In this work, we propose a posteriori error estimators and an adaptive reconstruction algorithm of the form (1.1) for (1.3) based on a separate marking using three error indicators in the module MARK; see algorithm 1. Further, we give a convergence analysis of the algorithm, in the sense that the sequence of state, adjoint and conductivity generated by the adaptive algorithm contains a convergent subsequence to a solution of the necessary optimality condition. The technical proof consists of two steps: step 1 shows the subsequential convergence to a solution of a limiting problem, and step 2 proves that the solution of the limiting problem satisfies the necessary optimality condition. The main technical challenges in the convergence analysis include the nonlinearity of the forward model, the nonconvexity of the double-well potential and properly treating the variational inequality. The latter two are overcome by pointwise convergence of discrete minimizers and Lebesgue's dominated convergence theorem, and AFEM analysis techniques for elliptic obstacle problems, respectively. The adaptive algorithm and its convergence analysis are the main contributions of this work.

Last, we situate this work in the existing literature. In recent years, several adaptive techniques, including AFEM, have been applied to the numerical resolution of inverse problems. In a series of works [6–9], Beilina et al studied the AFEM in a dual-weighted residual framework for parameter identification. Feng et al [20] proposed a residual-based estimator for the state, adjoint and control by assuming convexity of the cost functional and high regularity on the control. Li et al [38] derived a posteriori error estimators for recovering the flux and proved their reliability; see [55] for a plain convergence analysis. Clason et al [17] studied functional a posteriori estimators for convex regularized formulations. Recently, Jin et al [32] proposed a first AFEM for the Tikhonov functional for EIT with an $\newcommand{\Om}{\Omega} H^1(\Omega)$ penalty, and also provided a convergence analysis. This work extends the approach in [32] to the case of a piecewise constant conductivity. There are a number of major differences between this work and [32]. First, the $\newcommand{\Om}{\Omega} H^1(\Omega)$ penalty in [32] facilitates deriving the a posteriori estimator on the conductivity $\sigma$, by completing the squares and suitable approximation argument, which is not directly available for the Mordica–Mortola functional $\mathcal{F}_\varepsilon$. Second, we develop a sharper error indicator associated with the crucial variational inequality than that in [32], by a novel constraint-preserving interpolation operator [13]; see the proof of theorem 5.5, which represents the main technical novelty of this work. Third, algorithm 1 employs a separate marking for the estimators, instead of the collective marking in [32], which automatically takes care of different scalings of the estimators.

The rest of this paper is organized as follows. In section 2, we introduce the CEM and the regularized least-squares formulation. In section 3, we give the AFEM algorithm. In section 4, we present extensive numerical results to illustrate its convergence and efficiency. In section 5, we present the lengthy technical convergence analysis. Throughout, $\langle\cdot,\cdot\rangle$ and $(\cdot,\cdot)$ denote the inner product on the Euclidean space and $\newcommand{\Om}{\Omega} (L^2(\Omega)){}^d$, respectively, by $\|\cdot\|$ the Euclidean norm, and occasionally abuse $\langle\cdot,\cdot\rangle$ for the duality pairing between the Hilbert space $\mathbb{H}$ and its dual space. The superscript ${\rm t}$ denotes the transpose of a vector. The notation c denotes a generic constant, which may differ at each occurrence, but it is always independent of the mesh size and other quantities of interest.

This part describes the regularized approach for recovering piecewise constant conductivities.

2.1. Complete electrode model (CEM)

Let $\newcommand{\Om}{\Omega} \Omega$ be an open bounded domain in $\mathbb{R}^{d}$ $(d=2,3)$ with a polyhedral boundary $\newcommand{\Om}{\Omega} \partial\Omega$. We denote the set of electrodes by $\{e_l\}_{l=1}^L$, which are line segments/planar surfaces on $\newcommand{\Om}{\Omega} \partial\Omega$ and satisfy $\newcommand{\e}{{\rm e}} \bar{e}_i\cap\bar{e}_k=\emptyset$ if $i\neq k$. The applied current on electrode e_l is denoted by I_l, and the vector $I=(I_1,\ldots,I_L){}^\mathrm{t}\in\mathbb{R}^L$ satisfies $\sum\nolimits_{l=1}^LI_l=0$, i.e. $I\in \mathbb{R}_\diamond^L :=\{V\in \mathbb{R}^L: \sum\nolimits_{l=1}^LV_l=0\}$. The electrode voltage $U=(U_1,\ldots,U_L){}^\mathrm{t}$ is normalized, i.e. $U\in\mathbb{R}_\diamond^L$. Then the CEM reads: given the conductivity $\sigma$, positive contact impedances $\{z_l\}_{l=1}^L$ and input current $I\in\mathbb{R}_\diamond^L$, find $\newcommand{\Om}{\Omega} (u,U)\in H^1(\Omega)\otimes\mathbb{R}_\diamond^L$ such that [14, 50]

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \label{eqn:cem} \left\{\begin{array}{@{}ll@{}} -\nabla\cdot(\sigma\nabla u)=0 & {{\rm in}}~\Omega, \nonumber \\ u+z_l\sigma\frac{\partial u}{\partial n}=U_l& {{\rm on}}~e_l, l=1,2,\ldots,L, \nonumber \\ \int_{e_l}\sigma\frac{\partial u}{\partial n}\mathrm{d}s =I_l& {{\rm for}}~ l=1,2,\ldots, L, \nonumber \\ \sigma\frac{\partial u}{\partial n}=0&{{\rm on}}~ \partial\Omega\backslash\cup_{l=1}^Le_l. \end{array}\right. \nonumber \end{align} \tag{ 2.1 }$$

The inverse problem is to recover the conductivity $\sigma$ from a noisy version $U^\delta$ of the electrode voltage $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}U(\sigma^{\dagger})$ (for the exact conductivity $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\sigma^{\dagger}$) corresponding to one or multiple input currents.

Below the conductivity $\sigma$ is assumed to be piecewise constant, i.e. in the admissible set

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \mathcal{A}:=\{\sigma\in L^\infty(\Om):~\sigma=c_0+(c_1-c_0)\chi_{\Om_1}\}, \nonumber \end{align*}$$

where the constants $c_1>c_0>0$ are known, $\newcommand{\Om}{\Omega} \overline{\Om}_1\subset\Om$ is an unknown open set with a Lipschitz boundary and $\newcommand{\Om}{\Omega} \chi_{\Om_1}$ denotes its characteristic function. We denote by $\mathbb{H}$ the space $\newcommand{\Om}{\Omega} H^1(\Omega)\otimes \mathbb{R}_\diamond^L$ with its norm given by

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \|(u,U)\|_{\mathbb{H}}^2 = \|u\|_{H^1(\Omega)}^2 + \|U\|^2. \nonumber \end{align*}$$

A convenient equivalent norm on the space $\mathbb{H}$ is given below.

Lemma 2.1. On the space $\mathbb{H}$, the norm $\|\cdot\|_\mathbb{H}$ is equivalent to the norm $\|\cdot\|_{\mathbb{H},*}$ defined by

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \|(u,U)\|_{\mathbb{H},*}^2 = \|\nabla u\|_{L^2(\Omega)}^2 + \sum_{l=1}^L\|u-U_l\|_{L^2(e_l)}^2. \nonumber \end{align*}$$

The weak formulation of the model (2.1) reads [50]: find $(u,U)\in \mathbb{H}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:cemweakform} a(\sigma,(u,U),(v,V)) = \langle I,V\rangle \quad \forall (v,V)\in \mathbb{H}, \nonumber \end{align} \tag{ 2.2 }$$

where the trilinear form $a(\sigma,(u,U),(v,V))$ on $\mathcal{A}\times\mathbb{H}\times\mathbb{H}$ is defined by

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle a(\sigma,(u,U),(v,V)) = (\sigma \nabla u ,\nabla v) +\sum_{l=1}^Lz_l^{-1}(u-U_l,v-V_l)_{L^2(e_l)}, \nonumber \end{align*}$$

where $(\cdot,\cdot)_{L^2(e_l)}$ denotes the $L^2(e_l)$ inner product. For any $\sigma\in\mathcal{A}$, $\{z_l\}_{l=1}^L$ and $I\in \Sigma_\diamond^L$, the existence and uniqueness of a solution $(u,U)\in\mathbb{H}$ to (2.2) follows from lemma 2.1 and the Lax–Milgram theorem.

2.2. Regularized reconstruction

For numerical reconstruction with a piecewise constant conductivity, the TV penalty is popular. The conductivity $\sigma$ is assumed to be in the space $\newcommand{\Om}{\Omega} \mathrm{BV}(\Om)$ of bounded variation [4, 19], i.e.

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \mathrm{BV} (\Om) = \left\{v \in L^1(\Om): |v|_{\mathrm{TV}(\Om)}<\infty \right\}, \nonumber \end{align*}$$

equipped with the norm $\newcommand{\Om}{\Omega} \|v\|_{\mathrm{BV}(\Om)}=\|v\|_{L^1(\Om)}+|v|_{\mathrm{TV}(\Om)}$, where

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\dx}{\,{\rm d}x} \displaystyle |v|_{\mathrm{TV}(\Om)}:=\sup\left\{\int_\Om v \nabla\cdot\boldsymbol{\phi}\dx:~\boldsymbol{\phi}\in (C_c^1(\Om))^d,~\|\boldsymbol{\phi}(x)\|\leqslant 1\right\}. \nonumber \end{align*}$$

Below we discuss only one dataset, since the case of multiple datasets is similar. Then Tikhonov regularization leads to the following minimization problem:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\J}{\mathcal{J}} \newcommand{\Om}{\Omega} \displaystyle \label{eqn:tikh} \min_{\sigma\in\mathcal{A}}\left\{\J(\sigma) = \frac{1}{2} \|U(\sigma)-U^\delta\|^2 + \alpha|\sigma|_{{\rm TV}(\Om)}\right\}. \nonumber \end{align} \tag{ 2.3 }$$

The scalar $\alpha>0$ is known as a regularization parameter. It has at least one minimizer [22, 46].

Since $\sigma$ is piecewise constant, by Lebesgue decomposition theorem [4], the TV term $\newcommand{\Om}{\Omega} |\sigma|_{{\rm TV}(\Omega)}$ in (2.3) reduces to $\newcommand{\dd}{\,{\rm d}} \int_{S_{\sigma}} |[\sigma]| \dd \mathcal{H}^{d-1}$, where $S_\sigma$ is the jump set, $[\sigma]=\sigma^+-\sigma^-$ denotes the jump across $S_\sigma$ and $\mathcal{H}^{d-1}$ refers to the $(d-1)$-dimensional Hausdorff measure. The numerical approximation of (2.3) requires simultaneously treating two sets of different Hausdorff dimensions (i.e. $\newcommand{\Om}{\Omega} \Om$ and $S_\sigma$), which is very challenging. Thus, we replace the TV term $\newcommand{\Om}{\Omega} |\sigma|_{\rm TV(\Omega)}$ in (2.3) by a Modica–Mortola-type functional [43]:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\dx}{\,{\rm d}x} \displaystyle \mathcal{F}_\varepsilon(z):= \left\{\begin{array}{@{}ll@{}} \varepsilon \|\nabla z\|^2_{L^2(\Om)}+\frac{1}{\varepsilon}\int_{\Om}W(z)\dx & {\rm if}~z\in H^1(\Om), \nonumber \\ +\infty & {\rm otherwise}, \end{array} \right. \nonumber \end{align*}$$

where $\varepsilon$ is a small positive constant controlling the width of the transition interface, and $W: \mathbb{R}\rightarrow\mathbb{R}$ is the double-well potential given in (1.2). The functional $\mathcal{F}_\varepsilon$ was first proposed to model the phase transition of two immiscible fluids in [11]. It is connected with the TV semi-norm as follows [2, 42, 43]; see [10] for an introduction to $\Gamma$-convergence.

Theorem 2.1. With $\newcommand{\dd}{\,{\rm d}} c_W=\int_{c_0}^{c_1}\sqrt{W(s)}\dd s$, let

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \mathcal{F}(z):=\left\{\begin{array}{@{}ll@{}} 2c_W|z|_{\rm TV(\Om)} & {\rm if}~z\in \mathrm{BV}(\Om)\cap\mathcal{A}, \nonumber \\ +\infty & {\rm otherwise}. \end{array} \right. \nonumber \end{align*}$$

Then $\mathcal{F}_\varepsilon$ $\Gamma$ converges to $\mathcal{F}$ in $\newcommand{\Om}{\Omega} L^1(\Om)$ as $\varepsilon\rightarrow 0^+$. Let $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \{\eps_n\}_{n\geqslant1}$ and $\{v_{n}\}_{n\geqslant 1}$ be given sequences such that $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \eps_n\rightarrow 0^+$ and $\newcommand{\F}{{\mathcal F}} \newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \{\F_{\eps_n}(v_n)\}_{n\geqslant1}$ is bounded. Then $v_n$ is precompact in $\newcommand{\Om}{\Omega} L^1(\Om)$.

The proposed EIT reconstruction method reads

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\J}{\mathcal{J}} \displaystyle \label{eqn:tikh-MM} \inf_{\sigma\in \widetilde{\mathcal{A}}} \left\{\J_{\varepsilon}(\sigma) = \frac{1}{2} \|U(\sigma)-U^\delta\|^2 + \frac{\widetilde{\alpha}}{2}\mathcal{F}_\varepsilon(\sigma)\right\}, \nonumber \end{align} \tag{ 2.4 }$$

where $\widetilde{\alpha}=\alpha/c_W$, and the admissible set $\widetilde{\mathcal{A}}$ is defined as

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \widetilde{\mathcal{A}}:=\left\{{\sigma \in H^1(\Omega)}: c_0\leqslant \sigma(x)\leqslant c_1~ {\rm a.e.}~ x\in\Omega\right\}. \nonumber \end{align*}$$

Now we recall a useful continuity result of the forward map [22, lemma 2.2], which gives the continuity of the fidelity term in the functional $\mathcal{J}_\varepsilon$. See also [18, 30] for related continuity results.

Lemma 2.2. Let $\{\sigma_n\}_{n\geqslant1}\subset\widetilde{\mathcal{A}}$ satisfy $\sigma_n\rightarrow\sigma^\ast$ in $\newcommand{\Om}{\Omega} L^1(\Om)$. Then

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:fm_cont} \lim_{n\rightarrow\infty}\|\left(u(\sigma_n)-u(\sigma^\ast),U(\sigma_n)- U(\sigma^\ast)\right)\|_{\mathbb{H}}=0. \nonumber \end{align} \tag{ 2.5 }$$

Lemma 2.2 implies that $\{\mathcal{J}_\varepsilon\}_{\varepsilon>0}$ are continuous perturbations of $\mathcal{J}$ in $\newcommand{\Om}{\Omega} L^1(\Om)$. Then the stability of $\Gamma$-convergence [2, proposition 1(ii)] [10, remark 1.4] and theorem 2.1 indicate that $\newcommand{\J}{\mathcal{J}} \newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \J_\eps$ $\Gamma$ converges to $\mathcal{J}$ with respect to $\newcommand{\Om}{\Omega} L^1(\Om)$, and $\mathcal{J}_{\varepsilon}$ can (approximately) recover piecewise constant conductivities. Next we show the existence of a minimizer to $\newcommand{\J}{\mathcal{J}} \newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \J_\eps$.

Theorem 2.2. For each $\varepsilon>0$, there exists at least one minimizer to problem (2.4).

Proof. Since $\mathcal{J}_\varepsilon$ is non-negative, there exists a minimizing sequence $\{\sigma_n\}_{n\geqslant1}\subset \widetilde{\mathcal{A}}$ such that $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \mathcal{J}_\varepsilon(\sigma_n)\rightarrow m_\eps:=\inf\nolimits_{\sigma\in\widetilde{\mathcal{A}}}\mathcal{J}_\varepsilon(\sigma)$. Thus, $\newcommand{\Om}{\Omega} \sup\nolimits_n\|\nabla\sigma_n\|_{L^2(\Om)}<\infty$, which, along with $c_0\leqslant \sigma_n\leqslant c_1$, yields $\newcommand{\Om}{\Omega} \|\sigma_n\|_{H^1(\Om)}\leqslant c$. Since $\widetilde{\mathcal{A}}$ is closed and convex, there exists a subsequence, relabeled as $\{\sigma_n\}_{n\geqslant 1}$, and some $\sigma^\ast\in\widetilde{\mathcal{A}}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \label{pf:min_cont01} \sigma_{n}\rightharpoonup\sigma^\ast\quad {{\rm weakly~in}}~H^1(\Om),\quad \sigma_n\rightarrow\sigma^\ast\quad {{\rm in}}~L^1(\Om),\quad \sigma_n\rightarrow\sigma^\ast\quad {{\rm a.e.~in}}~\Om. \nonumber \end{align} \tag{ 2.6 }$$

Since $W(s)\in C^2[c_0,c_1]$, $\{W(\sigma_n)\}_{n\geqslant1}$ is uniformly bounded in $\newcommand{\Om}{\Omega} L^\infty(\Omega)$ and converges to $W(\sigma^\ast)$ almost everywhere in $\newcommand{\Om}{\Omega} \Om$. By Lebesgue's dominated convergence theorem [19, p 28, theorem 1.19], $\newcommand{\dx}{\,{\rm d}x} \newcommand{\Om}{\Omega} \int_\Om W(\sigma_n)\dx\rightarrow\int_\Om W(\sigma^\ast)\dx.$ By lemma 2.2 and the weak lower semi-continuity of the $\newcommand{\Om}{\Omega} H^1(\Omega)$ semi-norm, we obtain

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\eps}{\varepsilon} \newcommand{\J}{\mathcal{J}} \displaystyle \J_\varepsilon(\sigma^\ast)\leqslant \liminf_{n\rightarrow\infty}\J_{\varepsilon}(\sigma_n)\leqslant \lim_{n\rightarrow\infty}\J_{\varepsilon}(\sigma_n)=m_\eps. \nonumber \end{align*}$$

Thus $\sigma^\ast$ is a global minimizer of the functional $\mathcal{J}_\varepsilon$. □

To obtain the necessary optimality system of (2.4), we use the standard adjoint technique. The adjoint problem for (2.2) reads: find $(\,p,P)\in\mathbb{H}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:cem-adj} a(\sigma,(\,p,P),(v,V)) = \langle U(\sigma)-U^\delta,V\rangle \quad \forall (v,V)\in \mathbb{H}. \nonumber \end{align} \tag{ 2.7 }$$

By straightforward computation, the Gâteaux derivative $\mathcal{J}_\varepsilon'(\sigma)[\mu]$ of the functional $\mathcal{J}_{\varepsilon}$ at $\sigma\in\widetilde{\mathcal{A}}$ in the direction $\newcommand{\Om}{\Omega} \mu\in H^1(\Om)$ is given by

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\eps}{\varepsilon} \displaystyle \mathcal{J}_{\varepsilon}'(\sigma)[\mu] = \widetilde{\alpha}\left[\eps(\nabla\sigma,\nabla\mu) + \frac{1}{2\eps}(W'(\sigma),\mu)\right] - (\mu \nabla u(\sigma),\nabla p(\sigma)). \nonumber \end{align*}$$

Then the minimizer $\sigma^*$ to problem (2.4) and the respective state $(u^*,U^*)$ and adjoint $(\,p^*,P^*)$ satisfy the following necessary optimality system:

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:cem-optsys} \left\{\begin{array}{@{}l@{}} a(\sigma^*,(u^*,U^*),(v,V)) = \langle I,V\rangle \quad \forall (v,V)\in \mathbb{H}, \nonumber \\ a(\sigma^*,(\,p^*,P^*),(v,V)) = \langle U^*-U^\delta,V\rangle \quad \forall(v,V)\in\mathbb{H}, \nonumber \\ \widetilde{\alpha}\varepsilon (\nabla\sigma^{\ast},\nabla(\mu-\sigma^{\ast})) + \frac{\widetilde{\alpha}}{2\varepsilon}(W'(\sigma^\ast),\mu-\sigma^\ast) - ((\mu-\sigma^{\ast})\nabla u^{\ast},\nabla p^{\ast})\geqslant 0\quad \forall \mu\in\widetilde{\mathcal{A}}, \end{array}\right. \nonumber \end{align} \tag{ 2.8 }$$

where the variational inequality in the last line is due to the box constraint in the admissible set $\widetilde{\mathcal{A}}$. The optimality system (2.8) forms the basis of the adaptive algorithm and its convergence analysis.

Now we develop an AFEM for problem (2.4). Let $\newcommand{\cT}{{\mathcal T}} \cT_0$ be a shape-regular triangulation of $\newcommand{\Om}{\Omega} \overline{\Omega}$ into simplicial elements, each intersecting at most one electrode surface e_l, and $\mathbb{T}$ be the set of all possible conforming triangulations of $\newcommand{\Om}{\Omega} \overline{\Omega}$ obtained from $\newcommand{\cT}{{\mathcal T}} \cT_0$ by the successive use of bisection. Then $\mathbb{T}$ is uniformly shape-regular, i.e. the shape-regularity of any mesh $\mathcal{T}\in\mathbb{T}$ is bounded by a constant depending only on $\newcommand{\cT}{{\mathcal T}} \cT_0$ [44, 52]. Over any $\newcommand{\cT}{{\mathcal T}} \cT\in\mathbb{T}$, we define a continuous piecewise linear finite element space

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\cT}{{\mathcal T}} \displaystyle V_\cT = \left\{v\in C(\overline{\Omega}): v|_T\in P_1(T)\ \forall T\in\cT\right\}, \nonumber \end{align*}$$

where P₁(T) consists of all linear functions on T. The space $\newcommand{\cT}{{\mathcal T}} V_\cT$ is used for approximating the state u and adjoint p , and the discrete admissible set $\widetilde{\mathcal{A}}_\mathcal{T}$ for the conductivity is given by

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cT}{{\mathcal T}} \displaystyle \widetilde{\mathcal{A}}_\cT:= V_\cT \cap {\widetilde{\mathcal{A}}}. \nonumber \end{align*}$$

Given $\newcommand{\cT}{{\mathcal T}} \newcommand{\Atilde}{\widetilde{{\mathcal A}}} \sigma_\cT\in \Atilde_\cT$, the discrete analogue of problem (2.2) is to find $\newcommand{\cT}{{\mathcal T}} \newcommand{\e}{{\rm e}} (u_\cT,U_\cT) \in \mathbb{H}_\cT \equiv V_\cT \otimes \mathbb{R}^L_\diamond$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cT}{{\mathcal T}} \displaystyle \label{eqn:dispoly} a(\sigma_\cT,(u_\cT,U_\cT),(v_\cT,V)) = \langle I, V\rangle \quad \forall (v_{\cT},V)\in \mathbb{H}_\cT. \nonumber \end{align} \tag{ 3.1 }$$

Then we approximate problem (2.4) by minimizing the following functional over $\newcommand{\cT}{{\mathcal T}} \newcommand{\Atilde}{\widetilde{{\mathcal A}}} \Atilde_\cT$:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\eps}{\varepsilon} \newcommand{\cT}{{\mathcal T}} \displaystyle \label{eqn:discopt} J_{\eps,\cT}(\sigma_{\cT}) = \frac{1}{2}\|U_{\cT}(\sigma_{\cT})-U^\delta\|^2 + \frac{\widetilde{\alpha}}{2}\mathcal{F}_\varepsilon(\sigma_\cT). \nonumber \end{align} \tag{ 3.2 }$$

Then, similar to theorem 2.2, there exists at least one minimizer $\newcommand{\cT}{{\mathcal T}} \sigma_\cT^*$ to (3.2), and the minimizer $\newcommand{\cT}{{\mathcal T}} \sigma_{\cT}^{\ast}$ and the related state $\newcommand{\cT}{{\mathcal T}} (u^*_\cT,U^*_\cT)\in \mathbb{H}_\cT$ and adjoint $\newcommand{\cT}{{\mathcal T}} (\,p^*_\cT,P^*_\cT)\in\mathbb{H}_{\cT}$ satisfy

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Atilde}{\widetilde{{\mathcal A}}} \newcommand{\cT}{{\mathcal T}} \displaystyle \label{eqn:cem-discoptsys} \left\{\begin{array}{@{}l@{}} a(\sigma_\cT^*,(u^*_\cT,U_\cT^*),(v,V)) = \langle I,V\rangle \quad \forall(v,V)\in \mathbb{H}_{\cT}, \nonumber \\ a(\sigma_\cT^*,(\,p^*_\cT,P^*_\cT),(v,V)) = \langle U^*_\cT-U^\delta,V\rangle\quad \forall(v,V)\in\mathbb{H}_{\cT}, \nonumber \\ \widetilde{\alpha}\varepsilon(\nabla\sigma_{\cT}^{\ast},\nabla(\mu-\sigma_{\cT}^{\ast})) + \frac{\widetilde{\alpha}}{2\varepsilon}(W'(\sigma_\cT^\ast),\mu-\sigma_\cT^\ast) -((\mu-\sigma_{\cT}^{\ast})\nabla u_{\cT}^{\ast},\nabla p_{\cT}^{\ast})\geqslant 0\quad \forall\mu\in\Atilde_{\cT}. \end{array}\right. \nonumber \end{align} \tag{ 3.3 }$$

Further, $\newcommand{\cT}{{\mathcal T}} (u_{\cT}^\ast,U_{\cT}^\ast)$ and $\newcommand{\cT}{{\mathcal T}} (\,p_{\cT}^\ast,P_{\cT}^\ast)$ depend continuously on the problem data, i.e.

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\cT}{{\mathcal T}} \displaystyle \label{stab-discpolyadj} \|(u_{\cT}^\ast,U_{\cT}^\ast)\|_{\mathbb{H}}+\|(\,p_{\cT}^\ast,P_{\cT}^\ast)\|_{\mathbb{H}}\leqslant c(\|I\|+\|U^\delta\|), \nonumber \end{align} \tag{ 3.4 }$$

where the constant c can be made independent of $\alpha$ and $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \eps$.

To describe the error estimators, we first recall some useful notation. The collection of all faces (respectively all interior faces) in $\newcommand{\cT}{{\mathcal T}} \cT\in\mathbb{T}$ is denoted by $\newcommand{\cT}{{\mathcal T}} \mathcal{F}_{\cT}$ (respectively $\newcommand{\cT}{{\mathcal T}} \mathcal{F}_{\cT}^i$) and its restriction to the electrode $\bar{e}_{l}$ and $\newcommand{\Om}{\Omega} \partial\Omega\backslash\cup_{l=1}^Le_l$ by $\newcommand{\cT}{{\mathcal T}} \mathcal{F}_{\cT}^{l}$ and $\newcommand{\cT}{{\mathcal T}} \mathcal{F}_{\cT}^{c}$, respectively. A face/edge F has a fixed normal unit vector $\boldsymbol{n}_{F}$ in $\newcommand{\Om}{\Omega} \overline{\Omega}$ with $\boldsymbol{n}_{F}=\boldsymbol{n}$ on $\newcommand{\Om}{\Omega} \partial\Omega$. The diameter of any $\newcommand{\cT}{{\mathcal T}} T\in\cT$ and $\newcommand{\cT}{{\mathcal T}} F\in\mathcal{F}_{\cT}$ is denoted by $h_{T}:=|T|^{1/d}$ and $h_{F}:=|F|^{1/(d-1)}$, respectively. For the solution $\newcommand{\cT}{{\mathcal T}} (\sigma^{\ast}_{\cT},(u _{\cT}^{\ast},U_{\cT}^{\ast}), (\,p_{\cT}^{\ast},P_{\cT}^{\ast}))$ to problem (3.3), we define two element residuals for each element $\newcommand{\cT}{{\mathcal T}} T\in\cT$ and two face residuals for each face $\newcommand{\cT}{{\mathcal T}} F\in\mathcal{F}_{\cT}$ by

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cT}{{\mathcal T}} \displaystyle R_{T,1}(\sigma_{\cT}^{\ast},u^{\ast}_{\cT}) & = \nabla\cdot(\sigma_{\cT}^{\ast}\nabla u^{\ast}_{\cT}),\nonumber \\ R_{T,2}(\sigma_{\cT}^\ast,u^{\ast}_{\cT},p^{\ast}_{\cT})& =\frac{\widetilde{\alpha}}{2\varepsilon}W'(\sigma_\cT^\ast)-\nabla u^{\ast}_{\cT}\cdot\nabla p^{\ast}_{\cT},\nonumber \\ J_{F,1}(\sigma_{\cT}^{\ast},u^{\ast}_{\cT},U_{\cT}^{\ast}) &= \left\{\begin{array}{@{}lll@{}} [\sigma_{\cT}^{\ast}\nabla u_{\cT}^{\ast}\cdot\boldsymbol{n}_{F}]\quad & {\rm for} ~F\in\mathcal{F}_{\cT}^i, \nonumber \\ \sigma_{\cT}^{\ast}\nabla u_{\cT}^{\ast}\cdot\boldsymbol{n}+(u_{\cT}^{\ast}-U_{\cT,l}^{\ast})/z_{l}\quad & {\rm for} ~ F\in\mathcal{F}_{\cT}^l, \nonumber \\ \sigma_{\cT}^{\ast}\nabla u_{\cT}^{\ast}\cdot\boldsymbol{n}\quad & {\rm for} ~F\in\mathcal{F}_{\cT}^{c}, \end{array}\right.\nonumber \\ J_{F,2}(\sigma^{\ast}_{\cT}) &= \left\{\begin{array}{@{}lll@{}} \widetilde{\alpha}\varepsilon[\nabla\sigma_{\cT}^{\ast}\cdot\boldsymbol{n}_{F}]\quad & {\rm for} ~F\in\mathcal{F}_{\cT}^i, \nonumber \\ \widetilde{\alpha}\varepsilon\nabla\sigma_{\cT}^{\ast}\cdot\boldsymbol{n}\quad & {\rm for} ~F\in \mathcal{F}_\cT^l \cup \mathcal{F}_\cT^c, \end{array}\right. \nonumber \end{align*}$$

where $[\cdot]$ denotes the jump across interior face F. Then for any element $\newcommand{\cT}{{\mathcal T}} T\in \cT$, we define the following three error estimators

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cT}{{\mathcal T}} \displaystyle \eta_{\cT,1}^{2}(\sigma_{\cT}^{\ast},u_{\cT}^{\ast},U_{\cT}^{\ast},T) & :=h_{T}^{2}\|R_{T,1}(\sigma^{\ast}_{\cT},u^{\ast}_{\cT})\|_{L^{2}(T)}^{2} +\sum_{F\subset\partial T}h_{F}\|J_{F,1}(\sigma^{\ast}_{\cT},u^{\ast}_{\cT},U^{\ast}_{\cT})\|_{L^{2}(F)}^{2},\nonumber \\ \eta_{\cT,2}^{2}(\sigma_{\cT}^{\ast},p_{\cT}^{\ast},P_{\cT}^{\ast},T) &:=h_{T}^{2}\|R_{T,1}(\sigma^{\ast}_{\cT},p^{\ast}_{\cT})\|_{L^{2}(T)}^{2} +\sum_{F\subset\partial T}h_{F}\|J_{F,1}(\sigma^{\ast}_{\cT},p^{\ast}_{\cT},P^{\ast}_{\cT})\|_{L^{2}(F)}^{2},\nonumber \\ \eta_{\cT,3}^{q}(\sigma_{\cT}^{\ast},u_{\cT}^{\ast},p_{\cT}^{\ast},T) &:=h_{T}^{q}\|R_{T,2}(\sigma_{\cT}^\ast,u_{\cT}^{\ast},p_{\cT}^{\ast})\|^{q}_{L^{q}(T)}+\sum_{F\subset\partial T}h_{F}\|J_{F,2}(\sigma_{\cT}^{\ast})\|^{q}_{L^{q}(F),} \nonumber \end{align*}$$

with $q=d/(d-1)$. The estimator $\newcommand{\cT}{{\mathcal T}} \newcommand{\e}{{\rm e}} \eta_{\cT,1}(\sigma_{\cT}^{\ast},u_{\cT}^{\ast},U_{\cT}^{\ast},T)$ is identical to the standard residual error indicator for the direct problem: find $(\tilde u,\tilde U)\in \mathbb{H}$ such that

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cT}{{\mathcal T}} \displaystyle a(\sigma_\cT^*,(\tilde u,\tilde U),(v,V)) = \langle I,V\rangle, \quad \forall (v,V)\in \mathbb{H}. \nonumber \end{align*}$$

It differs from the direct problem in (2.8) by replacing the conductivity $\sigma^*$ with $\newcommand{\cT}{{\mathcal T}} \sigma^*_\cT$ instead, and is a perturbation of the latter case. The perturbation is vanishingly small in the event of the conjectured (subsequential) convergence $\newcommand{\cT}{{\mathcal T}} \sigma^*_\cT\rightarrow \sigma^*$. The estimator $\newcommand{\cT}{{\mathcal T}} \newcommand{\e}{{\rm e}} \eta_{\cT,2}(\sigma_{\cT}^{\ast},p_{\cT}^{\ast}, P_{\cT}^{\ast},T)$ admits a similar interpretation. These two estimators are essentially identical to that for the $\newcommand{\Om}{\Omega} H^1 (\Omega)$ penalty in [32], and we refer to [32, section 3.3] for a detailed heuristic derivation. The estimator $\newcommand{\cT}{{\mathcal T}} \newcommand{\e}{{\rm e}} \eta_{\cT,3}(\sigma_{\cT}^{\ast},u_{\cT}^{\ast},p_{\cT}^{\ast},T)$ is related to the variational inequality in the necessary optimality condition (2.8), and roughly provides a quantitative measure of how well it is satisfied. The estimator (including the exponent q) is motivated by the convergence analysis; see the proof of theorem 5.5 and remark 5.2 below. It represents the main new ingredient for problem (2.4), and differs from that for the $\newcommand{\Om}{\Omega} H^1(\Omega)$ penalty in [32].

Remark 3.1. The estimator $\newcommand{\e}{{\rm e}} \eta_{k,3}$ improves that in [32], i.e.

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cT}{{\mathcal T}} \displaystyle \eta_{\cT,3}^2(\sigma_{\cT}^{\ast},u_{\cT}^{\ast},p_{\cT}^{\ast},T) :=h_{T}^4\|R_{T,2}(\sigma_{\cT}^\ast,u_{\cT}^{\ast},p_{\cT}^{\ast})\|^2_{L^2(T)}+\sum_{F\subset\partial T}h_{F}^2\|J_{F,2}(\sigma_{\cT}^{\ast})\|^2_{L^2(F)}, \nonumber \end{align*}$$

in terms of the exponents on h_T and h_F. This improvement is achieved by a novel constraint-preserving interpolation operator defined in (5.13) below.

Now we can formulate an adaptive algorithm for (2.4); see algorithm 1. Below we indicate the dependence on the mesh $\newcommand{\cT}{{\mathcal T}} \cT_k$ by the subscript k, e.g. $\mathcal{J}_{\varepsilon,k}$ for $\newcommand{\cT}{{\mathcal T}} \mathcal{J}_{\varepsilon,\cT_k}$.

Algorithm 1. AFEM for EIT with a piecewise constant conductivity.

1: Specify an initial mesh $\newcommand{\cT}{{\mathcal T}} \cT_{0}$, and set the maximum number K of refinements.

2: for $k=0:K-1$ do

3: (SOLVE) Solve problem (3.1) and (3.2) over $\newcommand{\cT}{{\mathcal T}} \cT_{k}$ for $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} (\sigma_k^{\ast},(u_{k}^{\ast},U_{k}^{\ast}))\in\Atilde_{k}\times\mathbb{H}_{k}$ and (3.3) for $(\,p_k^\ast,P_k^\ast)\in \mathbb{H}_k$.

4: (ESTIMATE) Compute error indicators $\newcommand{\e}{{\rm e}} \eta_{k,1}^{2}(\sigma_{k}^{\ast},u_{k}^{\ast},U_{k}^{\ast})$, $\newcommand{\e}{{\rm e}} \eta_{k,2}^{2}(\sigma_{k}^{\ast},p_{k}^{\ast},P_{k}^{\ast})$ and $\newcommand{\e}{{\rm e}} \eta_{k,3}^{q}(\sigma_{k}^{\ast},u_{k}^{\ast},p_{k}^{\ast})$.

5: (MARK) Mark three subsets $\newcommand{\cT}{{\mathcal T}} \mathcal{M}_{k}^i\subseteq\cT_{k}$ ($i=1,2,3$) such that each $\mathcal{M}_k^i$ contains at least one element

$\newcommand{\cT}{{\mathcal T}} \widetilde{T}_k^i\in\cT_{k}$ ($i=1,2,3$) with the largest error indicator:

$$\begin{align} \newcommand{\cT}{{\mathcal T}} \newcommand{\e}{{\rm e}} \eta_{k,i}(\widetilde{T}_k^i)=\max\nolimits_{T\in\cT_{k}}\eta_{k,i}\end{align} \tag{ 3.5 }$$

Then $\mathcal{M}_{k}:=\mathcal{M}_{k}^1\cup\mathcal{M}_{k}^2\cup\mathcal{M}_{k}^3$.

6: (REFINE) Refine each element T in $\mathcal{M}_{k}$ by bisection to get $\newcommand{\cT}{{\mathcal T}} \cT_{k+1}$.

7: Check the stopping criterion.

8: end for

9: Output $(\sigma_k^*,(u_k^*,U_k^*),(\,p_k^*,P_k^*))$.

The MARK module selects a collection of elements in the mesh $\mathcal{T}_k$. The condition (3.5) covers several commonly used marking strategies, e.g. maximum, equidistribution, modified equidistribution and Dörfler's strategy [49, pp 962]. Compared with a collective marking in AFEM in [32], algorithm 1 employs a separate marking to select more elements for refinement in each loop, which leads to fewer iterations of the adaptive process. The error estimators may also be used for coarsening, which is relevant if the recovered inclusions change dramatically during the iteration. However, the convergence analysis below does not carry over to coarsening, and it will not be further explored.

Last, we give the main theoretical result: for each fixed $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \eps>0$, the sequence of discrete solutions $\{\sigma_k^\ast, (u_k^\ast, U_k^\ast), (\,p_k^\ast, P_k^\ast)\}_{k\geqslant0}$ generated by algorithm 1 contains a subsequence converging in $\newcommand{\Om}{\Omega} H^1(\Om)\times\mathbb{H}\times\mathbb{H}$ to a solution of system (2.8). The proof is lengthy and technical, and thus deferred to section 5.

Theorem 3.1. The sequence of discrete solutions $\{\sigma_{k}^\ast,(u_{k}^\ast,U_{k}^\ast),(\,p_{k}^\ast,P_{k}^\ast)\}_{k\geqslant0}$ by algorithm 1 contains a subsequence $\{\sigma_{k_j}^\ast,(u_{k_j}^\ast,U_{k_j}^\ast),(\,p_{k_j}^\ast, P_{k_j}^\ast)\}_{j\geqslant0}$ convergent to a solution $(\sigma^\ast,(u^\ast,U^\ast),(\,p^\ast,P^\ast))$ of system (2.8):

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \|\sigma^\ast_{k_j}-\sigma^\ast\|_{H^1(\Omega)},~\|(u^\ast_{k_j}-u^\ast,U^\ast_{k_j}-U^\ast)\|_{\mathbb{H}},~ \|(\,p^\ast_{k_j}-p^\ast,P^\ast_{k_j}-P^\ast)\|_{\mathbb{H}}\rightarrow 0\quad {\rm as}~j\rightarrow\infty. \nonumber \end{align*}$$

Now we present numerical results to illustrate algorithm 1 on a square domain $\newcommand{\Om}{\Omega} \Omega=(-1,1){}^2$. There are sixteen electrodes $\{e_l\}_{l=1}^L$ (with L  =  16) evenly distributed along $\newcommand{\Om}{\Omega} \partial\Omega$, each of length 1/4. The contact impedances $\{z_l\}_{l=1}^L$ are all set to unit. We take ten sinusoidal input currents, and for each voltage $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}U(\sigma^{\dagger})\in\mathbb{R}^L_\diamond$, generate the noisy data $U^\delta$ by

$$\begin{align} \renewcommand{\dag}{\dagger} \newcommand{\e}{{\rm e}} \newcommand{\re}{{\rm Re}} \newcommand{\eps}{\varepsilon} \displaystyle \label{eqn:noisydata} U^\delta_l = U_l(\sigma^{\dagger}) + \epsilon \max_l|U_l(\sigma^{\dagger})|\xi_l,\ \ l=1,\ldots, L, \nonumber \end{align} \tag{ 4.1 }$$

where $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \epsilon$ is the (relative) noise level, and $\{\xi_l\}_{l=1}^L$ follow the standard normal distribution. Note that $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \epsilon=1\times 10^{-2}$ refers to a relatively high noise level for EIT. The exact data $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}U(\sigma^{\dagger})$ is computed using a much finer uniform mesh, to avoid the most obvious form of 'inverse crime'.

In the experiments, we fix K (the number of refinements) at 15, q (exponent in $\newcommand{\e}{{\rm e}} \eta_{k,3}^q$) at 2, and $\varepsilon$ (the functional $\mathcal{F}_\varepsilon$) at $1\times 10^{-2}$. The marking strategy (3.5) in the module MARK selects a minimal refinement set $\newcommand{\cT}{{\mathcal T}} \mathcal{M}_k{:=\cup_{i=1}^3\mathcal{M}_{k}^i}\subseteq\cT_k$ such that

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \eta_{k,1}^2(\sigma_k^{\ast},u_k^{\ast},U_k^{\ast},\mathcal{M}_{k}^1) &\geqslant\theta \eta_{k,1}^2(\sigma_k^{\ast},u_k^{\ast},U_k^{\ast}),\quad \eta_{k,2}^2(\sigma_k^{\ast},p_k^{\ast},P_k^{\ast},\mathcal{M}_{k}^2)\geqslant\theta \eta_{k,2}^2(\sigma_k^{\ast},p_k^{\ast},P_k^{\ast}),\nonumber \\ &\eta_{k,3}^2(\sigma_k^{\ast},u_k^{\ast},p_k^{\ast},\mathcal{M}_{k}^3)\geqslant\theta \eta_{k,3}^2(\sigma_k^{\ast},u_k^{\ast},p_k^{\ast}), \nonumber \end{align*}$$

with a threshold $\theta=0.7$. The refinement is performed with one popular refinement strategy, i.e. newest vertex bisection [41]. Specifically, it connects the midpoint x_T, as a newest vertex, of a reference edge F of an element $\newcommand{\cT}{{\mathcal T}} T\in\cT_k$ to the opposite node of F, and employs two edges opposite to the midpoint x_T as reference edges of the two newly created triangles in $\newcommand{\cT}{{\mathcal T}} \cT_{k+1}$. Problems (3.1) and (3.2) are solved by a Newton-type method; see appendix A for the details. The conductivity on $\newcommand{\cT}{{\mathcal T}} \cT_0$ is initialized to $\sigma_0=c_0$, and then for $k=1,2,...$, $\sigma_{k-1}^*$ (defined on $\newcommand{\cT}{{\mathcal T}} \cT_{k-1}$) is interpolated to $\newcommand{\cT}{{\mathcal T}} \cT_k$ to warm-start the optimization. The regularization parameter $\widetilde\alpha$ in (2.4) is determined in a trial-and-error manner. All computations are performed using MATLAB 2018a on a personal laptop with 8.00 GB RAM and 2.5 GHz CPU.

The first set of examples is concerned with two inclusions.

Example 1. The background conductivity $\sigma_0(x)=1$.

• (i)  
  The true conductivity $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\sigma^{\dagger}$ is given by $\sigma_0(x)+\chi_{B_1}(x) + \chi_{B_2}(x)$, where B₁ and B₂ denote two circles centered at $(0,0.5)$ and $(0,-0.5)$, respectively, both with a radius 0.3.
• (ii)  
  The true conductivity $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\sigma^{\dagger}$ is given by $\sigma_0(x)+1 + 1.2{\rm e}^{-\frac{25(x_1^2+(x_2-0.5){}^2)}{2}} + 1.2{\rm e}^{-\frac{25(x_1^2+(x_2+0.5){}^2)}{2}}$, i.e. two Gaussian bumps centered at $(0,0.5)$ and $(0,-0.5)$.
• (iii)  
  The true conductivity $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\sigma^{\dagger}$ is given by $\sigma_0(x)+5\chi_{B_1}(x) + 5\chi_{B_2}(x)$, where B₁ and B₂ denote two circles centered at $(0,0.5)$ and $(0,-0.5)$, respectively, both with a radius of 0.3.

The numerical results for example 1(i) with $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \epsilon=1\times 10^{-3}$ and $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \epsilon=1\times 10^{-2}$ are shown in figures 1–5, where d.o.f. denotes the degree of freedom of the mesh. It is observed from figure 1 that with both uniform and adaptive refinements, the final recoveries have comparable accuracy and capture the inclusion locations well.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Next we examine the adaptive refinement process more closely. In figures 2 and 3, we show the meshes $\mathcal{T}_k$ during the iteration and the corresponding recoveries $\sigma_k$ for example 1(i) at two noise levels $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \epsilon = 1\times 10^{-3}$ and $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \epsilon=1\times 10^{-2}$, respectively. On the coarse mesh $\newcommand{\cT}{{\mathcal T}} \cT_0$, the recovery has very large errors and can only identify one component and thus fails to correctly identify the number of inclusions, due to the severe under-resolution of both state and conductivity. Nonetheless, algorithm 1 can correctly recover the two components with reasonable accuracy after several adaptive loops, and accordingly, the support of the recovery is gradually refined with its accuracy improving steadily. In particular, the inclusion locations stabilize after several loops, and thus coarsening of the mesh seems unnecessary. Throughout, the refinement occurs mainly in the regions around the electrode edges and internal interface, which is clearly observed for both noise levels. This is attributed to the separable marking strategy, which allows detecting different sources of singularities simultaneously. In figure 4, we display the evolution of the error indicators for example 1(i) with $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \epsilon=1\times 10^{-3}$. The estimators play different roles: $\newcommand{\e}{{\rm e}} \eta_{k,1}^2$ and $\newcommand{\e}{{\rm e}} \eta_{k,2}^2$ indicate the electrode edges during first iterations and then also the internal interface, whereas throughout $\newcommand{\e}{{\rm e}} \eta_{k,3}^2$ concentrates on the internal interface. Thus, $\newcommand{\e}{{\rm e}} \eta_{k,1}^2$ and $\newcommand{\e}{{\rm e}} \eta_{k,2}^2$ are most effective for resolving the state and adjoint, whereas $\newcommand{\e}{{\rm e}} \eta_{k,3}^2$ is effective for detecting internal jumps of the conductivity. The magnitude of $\newcommand{\e}{{\rm e}} \eta_{k,2}^2$ is much smaller than $\newcommand{\e}{{\rm e}} \eta_{k,1}^2$, since the boundary data $U^\delta-U(\sigma_k)$ for the adjoint are much smaller than the input current I for the state. Thus, a simple collective marking strategy (i.e. $\newcommand{\e}{{\rm e}} \eta_k^2 =\eta_{k,1}^2 +\eta_{k,2}^2 + \eta_{k,3}^2$) may miss the correct singularity, due to their drastically different scalings. In contrast, the separate marking in (3.5) can take care of the scaling automatically.

In figure 5, we plot the $\newcommand{\Om}{\Omega} L^2(\Omega)$ and $\newcommand{\Om}{\Omega} L^1(\Omega)$ errors of the recoveries versus d.o.f. N, where the recovery on the corresponding finest mesh is taken as the reference (since the recoveries by the adaptive and uniform refinements are slightly different; see figure 1). Due to the discontinuity of the sought-after conductivity, the $\newcommand{\Om}{\Omega} L^1(\Omega)$ norm is especially suitable for measuring the convergence. The convergence of the algorithm is clearly observed for both adaptive and uniform refinements. Further, with a fixed d.o.f., AFEM gives more accurate results than the uniform one in both error metrics. These observations show the computational efficiency of the adaptive algorithm.

Examples 1(ii) and (iii) are variations of example 1(i), and the results are presented in figures 6–9. The proposed approach assumes a piecewise constant conductivity with known lower and upper bounds. Example 1(ii) does not fulfill the assumption, since the true conductivity $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\sigma^{\dagger}$ is not piecewise constant. Thus the algorithm can only produce a piecewise constant approximation to the exact one. Nonetheless, the inclusion support is reasonably identified. When the noise level $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \epsilon$ increases from $1\times 10^{-3}$ to $1\times 10^{-2}$, the reconstruction accuracy deteriorates significantly; see figure 6. Example 1(iii) involves high-contrast inclusions, which are well known to be numerically more challenging. This is clearly observed in figure 8, where the recovery accuracy is inferior, especially for the noise level $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \epsilon=1\times 10^{-2}$. However, the adaptive refinement procedure works well similarly to the preceding examples: the refinement occurs mainly around the electrode edges and inclusion interface; see figures 7 and 9 for the details.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Now we consider one more challenging example with four inclusions.

Example 2. The true conductivity $\newcommand{\re}{{\rm Re}} \renewcommand{\dag}{\dagger}\sigma^{\dagger}$ is given by $\sigma_0(x)+\sum\nolimits_{i=1}^4\chi_{B_i}(x)$, with the circles B_i centered at $(0.6,\pm0.6)$ and $(-0.6,\pm0.6)$, respectively, all with a radius 0.2, and the background conductivity $\sigma_0(x)=1$.

The numerical results for example 2 are given in figures 10–12. The results are in excellent agreement with the observations from example 1: the algorithm converges steadily as the adaptive iteration proceeds, and with a low noise level, it can accurately recover all four inclusions, showing clearly the efficiency of the adaptive approach. The refinement is mainly around the electrode edge and interval interface.

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

Zoom In Zoom Out Reset image size

The lengthy and technical proof is divided into two steps: step 1 shows the convergence to an auxiliary minimization problem over a limiting admissible set in section 5.1, and step 2 shows that the solution of the auxiliary problem satisfies the necessary optimality system (2.8) in section 5.2. The overall proof strategy is similar to [32], and hence we omit the relevant arguments whenever possible.

5.1. Auxiliary convergence

Since the two sequences $\{\mathbb{H}_k\}_{k\geqslant0}$ and $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \{\Atilde_k\}_{k\geqslant0}$ generated by algorithm 1 are nested, we may define

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\Atilde}{\widetilde{{\mathcal A}}} \newcommand{\Om}{\Omega} \displaystyle \mathbb{H}_{\infty}:=\overline{\bigcup_{k\geqslant 0}\mathbb{H}_{k}}~ ({{\rm in}}~\mathbb{H}{\rm \mbox{-}norm})\quad {\rm and}\quad \Atilde_{\infty}:=\overline{\bigcup_{k\geqslant 0}\Atilde_{k}}~({\rm in}~H^{1}(\Omega){\rm \mbox{-}norm}). \nonumber \end{align*}$$

Clearly $\mathbb{H}_{\infty}$ is a closed subspace of $\mathbb{H}$. For the set $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \Atilde_{\infty}$, we have the following result [32, lemma 4.1].

Lemma 5.1. $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \Atilde_{\infty}$ is a closed convex subset of $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \Atilde$.

Over the limiting set $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \Atilde_{\infty}$, we define an auxiliary limiting minimization problem:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Atilde}{\widetilde{{\mathcal A}}} \newcommand{\J}{\mathcal{J}} \displaystyle \label{eqn:medopt} \min_{\sigma_{\infty}\in\Atilde_{\infty}}\left\{\J_{\varepsilon,\infty}(\sigma_{\infty}) = \frac{1}{2}\|U_{\infty}(\sigma_{\infty})-U^\delta\|^2 + \frac{\widetilde{\alpha}}{2}\mathcal{F}_\varepsilon(\sigma_\infty)\right\}, \nonumber \end{align} \tag{ 5.1 }$$

where $(u_{\infty},U_{\infty})\in \mathbb{H}_{\infty}$ satisfies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:medpoly} a(\sigma_\infty,(u_\infty,U_\infty),(v,V)) = \langle I,V\rangle \quad \forall (v,V)\in\mathbb{H}_\infty. \nonumber \end{align} \tag{ 5.2 }$$

By lemma 2.1 and the Lax–Milgram theorem, problem (5.2) is well posed for any fixed $\sigma_\infty\in\mathcal{A}_\infty$. The next result gives the existence of a minimizer to (5.1) and (5.2).

Theorem 5.1. There exists at least one minimizer to problem (5.1)–(5.2).

Proof. Let $\{\sigma_k^\ast, (u_k^\ast, U_k^\ast)\}_{k\geqslant 0}$ be the sequence of discrete solutions given by algorithm 1. Since $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} c_1\in \Atilde_k$ for all k, by (3.4), $\newcommand{\J}{\mathcal{J}} \newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \J_{\eps,k} (\sigma_k^\ast)\leqslant \J_{\eps,k}(c_1)\leqslant c$, and thus $\{\sigma_{k}^\ast\}_{k\geqslant 0}$ is uniformly bounded in $\newcommand{\Om}{\Omega} H^1(\Omega)$. By lemma 5.1 and Sobolev embedding, there exists a subsequence, denoted by $\{\sigma_{k_j}^\ast\}_{j\geqslant 0}$, and some $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \sigma^\ast\in \Atilde_\infty$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \label{pf:medmin_01} \sigma_{k_j}^\ast \rightharpoonup \sigma^\ast\quad {{\rm weakly~in}}~H^1(\Om),\quad \sigma_{k_j}^\ast \rightarrow \sigma^\ast\quad {{\rm in}}~L^2(\Om),\quad \sigma_{k_j}^\ast \rightarrow \sigma^\ast\quad {{\rm a.e.~in}}~\Om. \nonumber \end{align} \tag{ 5.3 }$$

Next we introduce a discrete analogue of problem (5.2) with $\sigma_\infty=\sigma^\ast$: find $(u_{k_j},U_{k_j})\in \mathbb{H}_{k_j}$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pf:medmin_02} a(\sigma^\ast,(u_{k_j},U_{k_j}),(v,V)) = \langle I,V\rangle \quad \forall (v,V)\in\mathbb{H}_{k_j}. \nonumber \end{align} \tag{ 5.4 }$$

By lemma 2.1, Cea's lemma and the construction of the space $\mathbb{H}_{\infty}$, the solution $(u^\ast_\infty,U_\infty^\ast)\in\mathbb{H}_\infty$ of (5.2) with $\sigma_\infty=\sigma^\ast$ satisfies

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pf:medmin_03} \|(u^\ast_\infty-u_{k_j},U^\ast_\infty-U_{k_j})\|_{\mathbb{H}}\leqslant c \inf_{(v,V)\in \mathbb{H}_{k_j}}\|(u^\ast_\infty-v,U^\ast_\infty-V)\|_{\mathbb{H}} \rightarrow 0. \nonumber \end{align} \tag{ 5.5 }$$

Taking the test function $(v,V)=(u_{k_j}-u_{k_j}^\ast,U_{k_j}-U_{k_j}^\ast)\in\mathbb{H}_{k_j}$ in the first line of (3.3) and (5.4) and then applying the Cauchy–Schwarz inequality lead to

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle & a(\sigma_{k_j}^\ast,(u_{k_j}-u_{k_j}^\ast,U_{k_j}-U_{k_j}^\ast),(u_{k_j}-u_{k_j}^\ast,U_{k_j}-U_{k_j}^\ast))\nonumber \\ & =((\sigma_{k_j}^\ast-\sigma^\ast)\nabla(u_{k_j}-u^\ast_{\infty}),\nabla(u_{k_j}-u_{k_j}^\ast)) +((\sigma_{k_j}^\ast-\sigma^\ast)\nabla u_{\infty}^\ast,\nabla(u_{k_j}-u_{k_j}^\ast))\nonumber \\ &\leqslant (\|(\sigma_{k_j}^\ast-\sigma^\ast)\nabla(u_{k_j}-u^\ast_{\infty})\|_{L^2(\Omega)} + \|(\sigma_{k_j}^\ast-\sigma^\ast)\nabla u_{\infty}^\ast\|_{L^2(\Omega)}) \|\nabla(u_{k_j}-u_{k_j}^\ast)\|_{L^2(\Omega)}. \nonumber \end{align*}$$

In view of (5.5), pointwise convergence in (5.3) and Lebesgue's dominated convergence theorem,

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \|(\sigma_{k_j}^\ast-\sigma^\ast)\nabla(u_{k_j}-u^\ast_{\infty})\|_{L^2(\Omega)}\leqslant c_1 \|\nabla(u_{k_j}-u^\ast_{\infty})\|_{L^2(\Omega)} \rightarrow 0,\quad \|(\sigma_{k_j}^\ast-\sigma^\ast)\nabla u_{\infty}^\ast\|_{L^2(\Omega)} \rightarrow 0. \nonumber \end{align*}$$

This and lemma 2.1 imply $\|(u_{k_j}-u_{k_j}^\ast,U_{k_j}-U_{k_j}^\ast)\|_{\mathbb{H}} \rightarrow 0.$ Then, (5.5) and the triangle inequality imply

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pf:medmin_04} \|(u_{k_j}^\ast-u_{\infty}^\ast,U_{k_j}^\ast-U_{\infty}^\ast)\|_{\mathbb{H}} \rightarrow 0. \nonumber \end{align} \tag{ 5.6 }$$

Meanwhile, repeating the argument of theorem 2.2 gives

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\dx}{\,{\rm d}x} \displaystyle \label{pf:medmin_05} \int_\Om W(\sigma_{k_j}^\ast)\dx\rightarrow \int_\Om W(\sigma^\ast)\dx. \nonumber \end{align} \tag{ 5.7 }$$

Next we apply a density argument. For any $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \sigma_{\infty}\in \Atilde_\infty$, by the construction of the space $\mathbb{H}_\infty$, there exists a sequence $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \newcommand{\bi}{\boldsymbol} \{\sigma_k\}_{k\geqslant 0}\subset \bigcup\nolimits_{k\geqslant0}\Atilde_{k}$ such that $\sigma_k\rightarrow\sigma_\infty$ in $\newcommand{\Om}{\Omega} H^1(\Om)$. Repeating the preceding argument gives $\|U(\sigma_{k})-U^\delta\|^2\rightarrow\|U(\sigma_\infty)-U^\delta\|^2$ and $\newcommand{\dx}{\,{\rm d}x} \newcommand{\Om}{\Omega} \int_{\Om} W(\sigma_k)\dx\rightarrow \int_\Om W(\sigma_{\infty})\dx$. Now (5.6), the weak lower semicontinuity of the $\newcommand{\Om}{\Omega} H^1(\Omega)$-norm, (5.7) and the minimizing property of $\sigma_{k}^\ast$ to $\newcommand{\J}{\mathcal{J}} \newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \J_{\eps,k}$ over the set $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \Atilde_k$ imply

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Atilde}{\widetilde{{\mathcal A}}} \newcommand{\eps}{\varepsilon} \newcommand{\J}{\mathcal{J}} \displaystyle \label{eqn:lsc-med} \J_{\eps,\infty}(\sigma^\ast) &\leqslant \liminf_{j\rightarrow\infty}\J_{\eps,k_j}(\sigma^\ast_{k_j})\leqslant \limsup_{j\rightarrow\infty}\J_{\eps,k_j}(\sigma^\ast_{k_j})\nonumber \\ &\leqslant\limsup_{k\rightarrow\infty}\J_{\eps,k}(\sigma^\ast_{k})\leqslant \limsup_{k\rightarrow\infty}\J_{\eps,k}(\sigma_{k})=\J_{\eps,\infty}(\sigma_\infty)\quad \forall\sigma_\infty\in \Atilde_{\infty}. \nonumber \end{align} \tag{ 5.8 }$$

Since $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \sigma^\ast\in\Atilde_{\infty}$, $\sigma_\infty^\ast:=\sigma^\ast$ is a minimizer of $\newcommand{\J}{\mathcal{J}} \newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \J_{\eps,\infty}$ over $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \Atilde_{\infty}$. □

Further, we have the following auxiliary convergence.

Theorem 5.2. The sequence of discrete solutions $\{\sigma_{k}^{\ast},(u^{\ast}_{k},U_{k}^{\ast})\}_{k\geqslant0}$ to problem (3.2) contains a subsequence $\{\sigma_{k_{j}}^{\ast},(u_{k_{j}}^{\ast},U_{k_{j}}^\ast)\}_{j\geqslant 0}$ convergent to a minimizer $(\sigma_{\infty}^{\ast},(u_{\infty}^{\ast},U_{\infty}^\ast))$ to problem (5.1)–(5.2):

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \sigma_{k_{j}}^\ast\rightarrow\sigma_{\infty}^\ast\quad {\rm in}~H^1(\Omega), \quad \sigma_{k_{j}}^\ast\rightarrow\sigma_{\infty}^\ast\quad {\rm a.e.~in}~\Omega, \quad (u_{k_{j}}^\ast,U_{k_j}^*)\rightarrow (u_{\infty}^\ast,U_\infty^*)\quad {\rm in}~\mathbb{H}. \nonumber \end{align*}$$

Proof. The convergence of $(u_{k_j}^*,U_{k_j}^*)$ was already proved in theorem 5.1. Taking $\sigma_\infty=\sigma_\infty^\ast$ in (5.8) gives $\newcommand{\J}{\mathcal{J}} \newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \lim_{j\rightarrow\infty}\J_{\eps,k_j}(\sigma_{k_j}^\ast) =\J_{\eps,\infty}(\sigma^\ast_\infty)$. By (5.6) and (5.7), we have $\newcommand{\Om}{\Omega} \|\nabla\sigma_{k_j}^\ast \|^2_{L^2(\Om)}\rightarrow \|\nabla\sigma_{\infty}^\ast\|^2_{L^2(\Om)}$. Thus, the sequence $\{\sigma_{k_j}^\ast\}_{j\geqslant 0}$ converges to $\sigma_{\infty}^\ast$ in $\newcommand{\Om}{\Omega} H^1(\Om)$. □

Next we consider the convergence of the sequence $\{(\,p_k^\ast,P_k^\ast)\}_{k\geqslant 0}$. With a minimizer $(\sigma_\infty^\ast,(u^\ast_\infty, U^\ast_\infty))$ to problem (5.1), we define a limiting adjoint problem: find $(\,p^\ast_{\infty},P^\ast_{\infty})\in\mathbb{H}_\infty$ such that

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:medadj} a(\sigma_\infty^*,(\,p^*_\infty,P_\infty^*),(v,V)) = \langle U_\infty^*-U^\delta,V\rangle\quad \forall (v,V)\in\mathbb{H}_\infty. \nonumber \end{align} \tag{ 5.9 }$$

By lemma 2.1 and Lax–Milgram theorem, (5.9) is uniquely solvable. We have the following convergence result for $(\,p_\infty^\ast, P_\infty^\ast)$. The proof is identical to [32, theorem 4.5], and hence omitted.

Theorem 5.3. Under the condition of theorem 5.2, the subsequence of adjoint solutions $\{(\,p_{k_{j}}^\ast,P_{k_{j}}^\ast)\}_{j\geqslant 0}$ generated by algorithm 1 converges to the solution $(\,p_{\infty}^\ast, P_{\infty}^\ast)$ of problem (5.9):

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \lim_{j\rightarrow\infty}\|(\,p_{k_j}^\ast-p_{\infty}^\ast,P_{k_j}^\ast-P_{\infty}^\ast)\|_{\mathbb{H}}=0. \nonumber \end{align*}$$

5.2. Proof of theorem 3.1

Theorem 3.1 follows directly by combining theorems 5.2 and 5.3 in section 5.1 and theorems 5.4 and 5.5 below. The proof in this part relies on the marking condition (3.5). First, we show that the limit $(\sigma^\ast_{\infty},(u^\ast_{\infty},U^\ast_{\infty}),(\,p^\ast_{\infty},P^\ast_{\infty}))$ solves the variational equations in (2.8).

Theorem 5.4. The solutions $(\sigma_\infty^*,u_\infty^*,U_\infty^*)$ and $(\,p_\infty^*,P_\infty^*)$ to problems (5.1), (5.2) and (5.9) satisfy

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle a(\sigma^*_\infty,(u^\ast_\infty,U_\infty^\ast),(v,V)) &= \langle I, V\rangle\quad \forall (v,V)\in\mathbb{H},\nonumber \\ a(\sigma^*_\infty,(\,p^\ast_\infty,P_\infty^\ast),(v,V)) &= \langle U_\infty^*-U^\delta, V\rangle\quad \forall (v,V)\in\mathbb{H}. \nonumber \end{align*}$$

Proof. The proof is identical to [32, lemma 4.8], using theorems 5.2 and 5.3, and hence we only give a brief sketch. By [32, lemma 3.5], for each $\newcommand{\cT}{{\mathcal T}} T\in\cT_k$ with its face F (intersecting with e_l), there hold

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \eta_{k,1}^2(\sigma_{k}^\ast,u^\ast_{k},U^\ast_{k},T) &\leqslant c(\|\nabla u^\ast_{k}\|^2_{L^2(D_{T})}+h_{F}\|u^\ast_{k}-U^{\ast}_{k,l}\|^2_{L^2(F\cap e_{l})}),\nonumber \\ \eta_{k,2}^2(\sigma_{k}^\ast,p^\ast_{k},P^\ast_{k},T)&\leqslant c(\|\nabla p^\ast_{k}\|^2_{L^2(D_{T})}+h_{F}\|p_k^\ast-P^{\ast}_{k,l}\|^2_{L^2(F\cap e_{l})}), \nonumber \end{align*}$$

where the notation D_T is defined below. Then by the marking condition (3.5), [32, lemma 4.6] implies that for each convergent subsequence $\{\sigma^{\ast}_{k_j},(u^{\ast}_{k_j},U^{\ast}_{k_j}),(\,p^{\ast}_{k_j}, P^{\ast}_{k_j})\}_{j\geqslant0}$ from theorems 5.2 and 5.3, there hold

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \lim_{j\rightarrow\infty}\max_{T\in\mathcal{M}_{k_j}^1}\eta_{k_{j},1}(\sigma_{k_j}^*,u_{k_j}^*,U_{k_j}^*,T)=0 \quad {\rm and}\quad \lim_{j\rightarrow\infty}\max_{T\in\mathcal{M}_{k_j}^2}\eta_{k_{j},2}(\sigma_{k_j}^*,p_{k_j}^*,P_{k_j}^*,T)=0. \nonumber \end{align*}$$

Last, by [32, lemma 4.7] and theorems 5.2 and 5.3, the argument of [32, lemma 4.8] completes the proof. □

Remark 5.1. The argument of theorem 5.4 dates back to [49], and the main tools include the Galerkin orthogonality of the residual operator, the Lagrange and the Scott–Zhang interpolation operators [16, 48], the marking condition (3.5) and a density argument. Further, the error estimators $\newcommand{\e}{{\rm e}} \eta_{k,1}(\sigma_k^*,u_k^*,U_k^*)$ and $\newcommand{\e}{{\rm e}} \eta_{k,2}(\sigma_k^*,p_k^*,P_k^*)$ emerge in the proof and are then employed in the module ESTIMATE of algorithm 1.

Next we prove that the limit $(\sigma^\ast_{\infty},(u^\ast_{\infty},U^\ast_{\infty}), (\,p^\ast_{\infty}, P^\ast_{\infty}))$ satisfies the variational inequality in (2.8). The proof relies crucially on a constraint-preserving interpolation operator. We denote by D_T the union of elements in $\newcommand{\cT}{{\mathcal T}} \cT$ with a non-empty intersection with an element $\newcommand{\cT}{{\mathcal T}} T\in\cT$, and by $\omega_F$ the union of elements in $\newcommand{\cT}{{\mathcal T}} \cT$ sharing a common face/edge with $\newcommand{\cT}{{\mathcal T}} F\in\mathcal{F}_\cT$. Let

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\cT}{{\mathcal T}} \displaystyle \cT_{k}^{+}:=\bigcap_{l\geqslant k}\cT_{l},\quad \cT_{k}^{0}:=\cT_{k}\setminus\cT_{k}^{+},\quad \Omega_{k}^{+}:=\bigcup_{T\in\cT^{+}_{k}}D_{T},\quad \Omega_{k}^{0}:=\bigcup_{T\in\cT^{0}_{k}}D_{T}. \nonumber \end{align*}$$

The set $\newcommand{\cT}{{\mathcal T}} \cT_{k}^{+}$ consists of all elements not refined after the kth iteration, and all elements in $\newcommand{\cT}{{\mathcal T}} \cT_{k}^{0}$ are refined at least once after the kth iteration. Clearly, $\newcommand{\cT}{{\mathcal T}} \cT_{l}^{+} \subset\cT_{k}^{+}$ for l  <  k. We also define a mesh-size function $\newcommand{\Om}{\Omega} h_{k}:\overline{\Omega}\rightarrow \mathbb{R}^{+}$ almost everywhere

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle h_k(x) = \left\{\begin{array}{@{}ll@{}} h_T, & \quad x\in T^i, \nonumber \\ h_F , & \quad x\in F^i, \end{array}\right. \nonumber \end{align*}$$

where Tⁱ denotes the interior of an element $\newcommand{\cT}{{\mathcal T}} T\in\cT_{k}$, and Fⁱ the relative interior of an edge $F\in\mathcal{F}_{k}$. It has the following property [49, corollary 3.3]:

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \label{eqn:conv_zero_mesh} \lim_{k\rightarrow\infty}\|h_{k}\chi_{\Omega_k^0}\|_{L^\infty(\Omega)}=0. \nonumber \end{align} \tag{ 5.10 }$$

The next result gives the limiting behavior of the maximal error indicator $\newcommand{\e}{{\rm e}} \eta_{k,3}$.

Lemma 5.2. Let $\{(\sigma^{\ast}_k,(u^{\ast}_k,U^{\ast}_k),(\,p^{\ast}_k,P^{\ast}_k))\}_{k\geqslant 0}$ be the sequence of discrete solutions generated by algorithm 1. Then for each convergent subsequence $\{\sigma^{\ast}_{k_j},(u^{\ast}_{k_j},U^{\ast}_{k_j}),(\,p^{\ast}_{k_j},P^{\ast}_{k_j})\}_{j\geqslant0}$, there holds

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \lim_{j\rightarrow\infty}\max_{T\in\mathcal{M}_{k_j}^3}\eta_{k_{j},3}(\sigma_{k_j}^*,u_{k_j}^*,p_{k_j}^*,T)=0. \nonumber \end{align*}$$

Proof. The inverse estimate and scaled trace theorem imply that for each $\newcommand{\cT}{{\mathcal T}} T\in\cT_k$ (with its face F)

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle h_{T}^q\|\frac{\tilde\alpha}{2\varepsilon}W'(\sigma_k^*) - \nabla u^*_k\cdot\nabla p^*_k\|^q_{L^q(T)} & \leqslant ch_T^q\|\nabla u^*_k\cdot\nabla p^*_k\|_{L^q(T)}^q+ ch_T^q\|W'(\sigma_k^*)\|_{L^q(T)}^q \nonumber \\ &\leqslant ch_T^{q}h_T^{d-dq}\|\nabla u^*_k\cdot\nabla p^*_k\|^q_{L^1(T)}+ch_T^q\|W'(\sigma_k^*)\|_{L^q(T)}^q,\nonumber \\ \sum_{F\subset\partial T}h_F\|J_{F,2}(\sigma_k^*)\|_{L^q(F)}^q&\leqslant c\sum_{F\subset\partial T}h_Fh_{F}^{-1}\|\nabla\sigma_{k}^*\|^q_{L^q(\omega_F)}. \nonumber \end{align*}$$

With the choice $q=d/(d-1)$, combining these two estimates gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:stab-indicator} \eta_{k,3}^q(\sigma_{k}^\ast,u^\ast_{k},p^\ast_{k},T)\leqslant c(\|\nabla u^*_k\cdot\nabla p^*_k\|^q_{L^1(T)}+h_T^q\|W'(\sigma_k^*)\|_{L^q(T)}^q+\|\nabla\sigma_{k}^\ast\|^q_{L^q(D_{T})}), \nonumber \end{align} \tag{ 5.11 }$$

where c depends on $\widetilde{\alpha}$ and $\newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \eps$ in $\newcommand{\F}{{\mathcal F}} \newcommand{\eps}{\varepsilon} \newcommand{\e}{{\rm e}} \F_\eps$. Next, for the subsequence $\{\sigma^{\ast}_{k_j},(u^{\ast}_{k_j},U^{\ast}_{k_j}),(\,p^{\ast}_{k_j},P^{\ast}_{k_j})\}_{j\geqslant0}$, let $\widetilde{T}_{j}^3 \in \mathcal{M}_{k_{j}}^3$ be the element with the largest error indicator $\newcommand{\e}{{\rm e}} \eta_{k_j,3}(\sigma_{k_j}^*,u_{k_j}^*,p_{k_j}^*,T)$. Since $\newcommand{\Om}{\Omega} D_{\widetilde{T}_j^i}\subset\Omega_{k_j}^0$, (5.10) implies

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \label{eqn:estmarked} |D_{\widetilde{T}^i_j}|\leqslant c\|h_{k_{j}}\|^d_{L^{\infty}(\Omega_{k_j}^0)}\rightarrow 0 \quad {{\rm as}}~j\rightarrow\infty. \nonumber \end{align} \tag{ 5.12 }$$

By (5.11), Cauchy–Schwarz inequality and triangle inequality, there holds

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \eta_{k_j,3}^q(\sigma^{\ast}_{k_j},&u^\ast_{k_j},p^\ast_{k_j},\widetilde{T}_{j}^3) \leqslant c(\|\nabla u^\ast_{k_j}\|_{L^2(\widetilde{T}_{j}^3)}^q\|\nabla p^\ast_{k_j}\|_{L^2(\widetilde{T}_{j}^3)}^q +h_{\widetilde{T}_{j}^3}^q\|W'(\sigma_k^*)\|_{L^q(\widetilde{T}_{j}^3)}^q+\|\nabla\sigma_{k_j}^\ast\|^q_{L^q(D_{\widetilde{T}_{j}^3})})\nonumber \\ &\leqslant c\big((\|\nabla (u^\ast_{k_j}-u^\ast_{\infty})\|_{L^2(\Om)}^q+\|\nabla u^\ast_{\infty}\|_{L^2(\widetilde{T}_{j}^3)}^q)(\|\nabla (\,p^\ast_{k_j}-p^\ast_{\infty})\|_{L^2(\Om)}^q+\|\nabla p^\ast_{\infty}\|_{L^2(\widetilde{T}_{j}^3)}^q)\nonumber \\ &\quad +h_{\widetilde{T}_j^3}^q(\|W'(\sigma_{k_j}^\ast)-W'(\sigma_\infty^\ast)\|^q_{L^q(\Om)}+\|W'(\sigma_\infty^\ast)\|^q_{L^q(\widetilde{T}_{j}^3)})\nonumber \\ &\quad +(\|\nabla(\sigma_{k_j}^\ast-\sigma_{\infty}^\ast)\|^q_{L^q(\Om)} +\|\nabla\sigma_{\infty}^\ast\|^q_{L^q(D_{\widetilde{T}_{j}^3})})\big). \nonumber \end{align*}$$

By theorems 5.2 and 5.3, Lebesgue's dominated convergence theorem, the choice $q=d/(d-1)\leqslant 2$ and Hölder inequality, we obtain $\newcommand{\Om}{\Omega} \|W'(\sigma_{k_j}^\ast)-W'(\sigma_\infty^\ast)\|^q_{L^q(\Om)}\rightarrow 0$ and $\newcommand{\Om}{\Omega} \|\nabla(\sigma_{k_j}^\ast- \sigma_{\infty}^\ast)\|^q_{L^q(\Om)}\rightarrow 0$. Then the absolute continuity of the norm $\newcommand{\Om}{\Omega} \|\cdot\|_{L^q(\Omega)}$ with respect to Lebesgue measure and (5.12) complete the proof. □

Due to a lack of Galerkin orthogonality for variational inequalities, we employ a local L^r-stable interpolation operator of Clément/Chen–Nochetto type. Let $\mathcal{N}_k$ be the set of all interior nodes of $\newcommand{\cT}{{\mathcal T}} \cT_k$ and $\{\phi_{x}\}_{x\in\mathcal{N}_k}$ be the nodal basis functions in $V_k$. For each $x\in\mathcal{N}_k$, the support of $\phi_x$ is denoted by $\omega_x$, i.e. the union of all elements in $\newcommand{\cT}{{\mathcal T}} \cT_k$ with a non-empty intersection with x. Then we define $\newcommand{\Om}{\Omega} \Pi_k:L^1(\Om)\rightarrow V_k$ by

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dx}{\,{\rm d}x} \displaystyle \label{eqn:cn_int_def} \Pi_k v := \sum_{x\in\mathcal{N}_k} \frac{1}{|\omega_x|}\int_{\omega_x}v \dx \phi_x. \nonumber \end{align} \tag{ 5.13 }$$

Clearly, $\newcommand{\Atilde}{\widetilde{{\mathcal A}}} \Pi_k v\in \Atilde_k$ if $c_0 \leqslant v \leqslant c_1$ a.e. $\newcommand{\Om}{\Omega} x\in\Om$. The definition is adapted from [13] (for elliptic obstacle problems) by replacing the maximal ball $\Delta_x\subset\omega_x$ centered at an interior node x by $\omega_x$. $\Pi_k v$ satisfies the following properties; see appendix B for a proof.

Lemma 5.3. For any $\newcommand{\Om}{\Omega} v\in W^{1,r}(\Omega)$, there hold for all $r\in[1, +\infty]$, any $\newcommand{\cT}{{\mathcal T}} T\in \cT_k$ and any $F\subset \partial T$,

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \|\Pi_k v\|_{L^r(T)} \leqslant c \|v\|_{L^r(D_T)}, \quad \|\nabla\Pi_k v\|_{L^r(T)} \leqslant c \|\nabla v\|_{L^r(D_T)},\nonumber \\ \| v - \Pi_k v \|_{L^r(T)} \leqslant c h_T \| \nabla v \|_{L^r(D_T)}, \quad \| v - \Pi_k v\|_{L^r(F)} \leqslant c h_F^{1-1/r} \|\nabla v\|_{L^r(D_T)}. \nonumber \end{align*}$$

Last we show that the limit $(\sigma^\ast_{\infty},(u^\ast_{\infty},U^\ast_{\infty}),(\,p^\ast_{\infty},P^\ast_{\infty}))$ satisfies the variational inequality in (2.8).

Theorem 5.5. The solutions $(\sigma_\infty^*,u_\infty^*,U_\infty^*)$ and $(\,p_\infty^*,P_\infty^*)$ to problems (5.1), (5.2) and (5.9) satisfy

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Atilde}{\widetilde{{\mathcal A}}} \newcommand{\eps}{\varepsilon} \displaystyle \widetilde{\alpha}\eps(\nabla\sigma^\ast_\infty,\nabla(\mu-\sigma^\ast_\infty)) + \frac{\widetilde{\alpha}}{2\eps}(W'(\sigma^\ast_\infty),\mu-\sigma^\ast_\infty) - (\nabla u_\infty^\ast,\nabla p_{\infty}^\ast(\mu-\sigma_{\infty}^\ast))\geqslant0\quad \forall\mu\in\Atilde. \nonumber \end{align*}$$

Proof. The proof is lengthy, and we break it into five steps.

Step i. Derive a preliminary variational inequality. We relabel the subsequence $\{\sigma_{k_j}^\ast,(u_{k_j}^\ast,U_{k_j}^\ast),(\,p_{k_j}^\ast,P_{k_j}^\ast)\}_{j\geqslant0}$ in theorems 5.2 and 5.3 as $\{\sigma_{k}^\ast,(u_{k}^\ast,U_{k}^\ast),(\,p_{k}^\ast,P_{k}^\ast)\}_{k\geqslant0}$. Let I_k be the Lagrange interpolation operator on $V_k$, and let $\alpha' = \widetilde\alpha\varepsilon$ and $\alpha''=\frac{\widetilde\alpha}{2\varepsilon}$. For any $\newcommand{\Om}{\Omega} \mu\in\widetilde{\mathcal{A}}\cap C^\infty(\overline{\Omega})$, $I_{k}\mu\in\widetilde{\mathcal{A}}_{k}$ and let $\nu=\mu-I_k\mu$. Direct computation gives

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \alpha'&(\nabla\sigma^\ast_k,\nabla(\mu-\sigma^\ast_k)) +\alpha''(W'(\sigma^\ast_k),\mu-\sigma^\ast_k) -((\mu-\sigma_{k}^\ast)\nabla u_k^\ast,\nabla p_{k}^\ast)\nonumber \\ &=\alpha'(\nabla\sigma^\ast_k,\nabla(\mu-I_k\mu)) +\alpha''(W'(\sigma^\ast_k),\mu-I_k\mu) -((\mu-I_{k}\mu)\nabla u_k^\ast,\nabla p_{k}^\ast)\nonumber \\ &\qquad +\alpha'(\nabla\sigma^\ast_k,\nabla(I_k\mu-\sigma_k^\ast)) +\alpha''(W'(\sigma^\ast_k),I_k\mu-\sigma^\ast_k) -((I_{k}\mu-\sigma_k^\ast)\nabla u_k^\ast,\nabla p_{k}^\ast)\nonumber \\ & =\alpha'(\nabla\sigma^\ast_k,\nabla(\nu-\Pi_k\nu)) +\alpha''(W'(\sigma^\ast_k),\nu-\Pi_k\nu) -((\nu-\Pi_{k}\nu)\nabla u_k^\ast,\nabla p_{k}^\ast) \nonumber \\ &\qquad +\alpha'(\nabla\sigma^\ast_k,\nabla\Pi_{k}\nu) +\alpha''(W'(\sigma^\ast_k),\Pi_{k}\nu) -(\Pi_k\nu\nabla u_k^\ast,\nabla p_{k}^\ast)\nonumber \\ &\qquad +\alpha'(\nabla\sigma^\ast_k,\nabla(I_k\mu-\sigma_k^\ast)) +\alpha''(W'(\sigma^\ast_k),I_k\mu-\sigma^\ast_k) -((I_{k}\mu-\sigma_k^\ast)\nabla u_k^\ast,\nabla p_{k}^\ast)\nonumber \\ &\geqslant\left[\alpha'(\nabla\sigma^\ast_k,\nabla(\nu-\Pi_k\nu)) +\alpha''(W'(\sigma^\ast_k),\nu-\Pi_k\nu) -((\nu-\Pi_{k}\nu)\nabla u_k^\ast,\nabla p_{k}^\ast)\right]\nonumber \\ &\quad +\left[\alpha'(\nabla\sigma^\ast_k,\nabla\Pi_{k}\nu) +\alpha''(W'(\sigma^\ast_k),\Pi_{k}\nu) -(\Pi_k\nu\nabla u_k^\ast,\nabla p_{k}^\ast)\right]:={\rm I}+ {\rm II},\label{eqn:gat_mc_01} \nonumber \end{align} \tag{ 5.14 }$$

where the last inequality is due to the variational inequality in (3.3) with $\mu_k=I_k\mu$.

Step ii. Bound the term ${\rm I}$. By elementwise integration by parts, Hölder inequality, the definition of the estimator $\newcommand{\e}{{\rm e}} \eta_{k,3}$ and lemma 5.3 with $r=q'$ (with $q'$ being the conjugate exponent of q),

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cT}{{\mathcal T}} \newcommand{\dx}{\,{\rm d}x} \displaystyle |{\rm I}| &= \left| \sum_{T \in \cT_k} \int_T R_{T,2}(\sigma_{k}^\ast, u_k^\ast, p_k^\ast)(\nu-\Pi_k\nu)\dx + \sum_{F\in\mathcal{F}_k}\int_F J_{F,2}(\sigma_k^\ast)(\nu-\Pi_k\nu) \mathrm{d}s \right| \nonumber \\ &\leqslant \sum_{T \in \cT_k}\Big(\|R_{T,2}(\sigma_{k}^\ast, u_k^\ast, p_k^\ast)\|_{L^q(T)}\|\nu-\Pi_k\nu\|_{L^{q'}(T)}+\sum_{F\subset\partial T}\| J_{F,2}(\sigma_k^\ast)\|_{L^{q}(F)}\|\nu-\Pi_k\nu\|_{L^{q'}(F)}\Big)\nonumber \\ &\leqslant c\sum_{T \in \cT_k}\Big(h_T\|R_{T,2}(\sigma_{k}^\ast, u_k^\ast, p_k^\ast)\|_{L^q(T)}+\sum_{F\subset\partial T}h_F^{1/q}\| J_{F,2}(\sigma_k^\ast)\|_{L^{q}(F)}\Big)\|\nabla\nu\|_{L^{q'}(D_T)}\nonumber \\ &\leqslant c\sum_{T\in\cT_{k}}\eta_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\|\nabla\nu\|_{L^{q'}(D_T)}. \nonumber \end{align*}$$

Thus, for any k  >  l, by (discrete) Hölder's inequality and the finite overlapping property of the patches D_T, due to uniform shape regularity of the meshes $\mathcal{T}_k\in\mathbb{T}$, there holds

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\cT}{{\mathcal T}} \displaystyle |{\rm I}| &\leqslant c\big(\sum_{T\in\cT_{k}\setminus\cT_l^+}\eta_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\|\nabla\nu\|_{L^{q'}(D_T)} +\sum_{T\in\cT_l^+}\eta_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\|\nabla\nu\|_{L^{q'}(D_T)}\big)\nonumber \\ &\leqslant c\Big(\big(\sum_{T\in\cT_{k}\setminus\cT^+_l}\eta^q_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\big)^{1/q}\|\nabla(\mu-I_k\mu)\|_{L^{q'}(\Omega_l^0)}\nonumber \\ &\qquad +\big(\sum_{T\in\cT_{l}^+}\eta^q_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\big)^{1/q}\|\nabla(\mu-I_k\mu)\|_{L^{q'}(\Omega_l^+)}\Big). \nonumber \end{align*}$$

Since $W'(s)\in C^1[c_0,c_1]$, by the pointwise convergence of $\{\sigma_{k}^\ast\}_{k\geqslant0}$ in theorem 5.2 and Lebesgue's dominated convergence theorem, we deduce

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \label{eqn:gat_mc_aux} W'(\sigma_{k}^\ast)\rightarrow W'(\sigma_\infty^\ast)\quad {{\rm in}}~L^2(\Om). \nonumber \end{align} \tag{ 5.15 }$$

Since $q=d/{(d-1)}\leqslant 2$, the sequence $\{W'(\sigma_k^\ast)\}_{k\geqslant0}$ is uniformly bounded in $\newcommand{\Om}{\Omega} L^q(\Omega)$. By theorems 5.2 and 5.3, the sequences $\{\sigma_{k}^\ast\}_{k\geqslant0}$, $\{u_{k}^\ast\}_{k\geqslant0}$ and $\{\,p_{k}^\ast\}_{k\geqslant0}$ are uniformly bounded in $\newcommand{\Om}{\Omega} H^1(\Omega)$. Thus, (5.11) and (5.10), and Hölder inequality give

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\cT}{{\mathcal T}} \displaystyle &\sum_{T\in\cT_{k}\setminus\cT^+_l}\eta^q_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\nonumber \\ &\leqslant c\Big(\|\nabla u_k^\ast\cdot\nabla p_k^\ast\|_{L^1(\Omega)}^{q-1}\sum_{T\in\cT_{k}\setminus\cT^+_l}\|\nabla u_k^\ast\cdot\nabla p_k^\ast\|_{L^1(T)} +\|h_{l}\|^{q}_{L^\infty(\Omega^{0}_l)}\|W'(\sigma_k^\ast)\|^q_{L^q(\Om)} +\|\nabla\sigma_{k}^\ast\|^q_{L^q(\Omega)}\Big)\nonumber \\ &\leqslant c\Big(\|\nabla u_{k}^\ast\|^q_{L^2(\Omega)}\|\nabla p_k^\ast\|^{q}_{L^2(\Omega)}+\|h_{l}\chi_{\Omega_l^0}\|^{q}_{L^\infty(\Omega)}\|W'(\sigma_k^\ast)\|^q_{L^q(\Om)} +\|\nabla\sigma_{k}^\ast\|^q_{L^2(\Omega)}\Big)\leqslant c.\label{eqn:bdd-eta} \nonumber \end{align} \tag{ 5.16 }$$

Then by the error estimate of I_k [16],

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\cT}{{\mathcal T}} \displaystyle |{\rm I}|\leqslant c\|h_{l}\chi_{\Omega^0_l}\|_{L^\infty(\Omega)}\|\mu\|_{W^{2,q'}(\Om)}+c\big(\sum_{T\in\cT_{l}^+}\eta^q_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\big)^{1/q}\|\mu\|_{W^{2,q'}(\Om)}. \nonumber \end{align*}$$

By (5.10), $\newcommand{\Om}{\Omega} c\|h_{l}\chi_{\Omega^0_l}\|_{L^\infty(\Omega)}\|\mu\|_{W^{2,q'}(\Om)}\rightarrow 0$ as $l\rightarrow\infty$. Since $\newcommand{\cT}{{\mathcal T}} \cT_l^+\subset\cT_{k}$ for k  >  l, (3.5) implies

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\cT}{{\mathcal T}} \displaystyle \big(\sum_{T\in\cT^+_l}\eta^q_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\big)^{1/q} \leqslant{|\cT_{l}^+|}^{1/q}\max_{T\in\cT_l^+}\eta_{k,3}(\sigma_{k}^\ast,u_{k}^\ast,p_{k}^\ast,T)\leqslant{|\cT^+_l|}^{1/q}\max_{T\in\mathcal{M}_k^3}\eta_{k,3}(\sigma_{k}^\ast,u_{k}^\ast,p_{k}^\ast,T). \nonumber \end{align*}$$

By lemma 5.2, for any small $\varepsilon>0$, we can choose $k_1>l_1$ for some large fixed l₁ such that whenever k  >  k₁,

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\cT}{{\mathcal T}} \displaystyle c({\sum_{T\in\cT_{l}^+}\eta^q_{k,3}}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T))^{1/q}\|\mu\|_{W^{2,q'}(\Omega)}<\varepsilon. \nonumber \end{align*}$$

Consequently,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \label{eqn:gat_mc_02} {\rm I} \rightarrow 0\quad \forall\mu\in\widetilde{\mathcal{A}}\cap C^\infty(\overline{\Om}). \nonumber \end{align} \tag{ 5.17 }$$

Step iii. Bound the term ${\rm II}$. For the term ${\rm II}$, elementwise integration and Hölder inequality yield

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\cT}{{\mathcal T}} \newcommand{\dx}{\,{\rm d}x} \displaystyle |{\rm II}|&= \left| \sum_{T \in \cT_k} \int_T R_{T,2}(\sigma_{k}^\ast, u_k^\ast, p_k^\ast)\Pi_k\nu\dx + \sum_{F\in\mathcal{F}_k}\int_F J_{F,2}(\sigma_k^\ast)\Pi_k\nu \mathrm{d}s \right|\nonumber \\ &\leqslant \sum_{T \in \cT_k}\Big(\|R_{T,2}(\sigma_{k}^\ast, u_k^\ast, p_k^\ast)\|_{L^q(T)}\|\Pi_k\nu\|_{L^{q'}(T)}+\sum_{F\subset\partial T}\| J_{F,2}(\sigma_k^\ast)\|_{L^{q}(F)}\|\Pi_k\nu\|_{L^{q'}(F)}\Big). \nonumber \end{align*}$$

By the scaled trace theorem, local inverse estimate, $L^{q'}$-stability of $\Pi_k$ in lemma 5.3, local quasi-uniformity and interpolation error estimate for I_k [16], we deduce that for k  >  l

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\cT}{{\mathcal T}} \displaystyle |{\rm II}|& \leqslant c\sum_{T \in \cT_k}\Big(h_T\|R_{T,2}(\sigma_{k}^\ast, u_k^\ast, p_k^\ast)\|_{L^q(T)}h_T^{-1}\|\Pi_k\nu\|_{L^{q'}(T)}+\sum_{F\subset\partial T}h_F^{1/q}\| J_{F,2}(\sigma_k^\ast)\|_{L^{q}(F)}h_F^{-1/q-1/q'}\|\Pi_k\nu\|_{L^{q'}(T)}\Big) \nonumber \\ &\leqslant c\sum_{T \in \cT_k}\Big(h_T\|R_{T,2}(\sigma_{k}^\ast, u_k^\ast, p_k^\ast)\|_{L^q(T)}+\sum_{F\subset\partial T}h_F^{1/q}\| J_{F,2}(\sigma_k^\ast)\|_{L^{q}(F)}\Big)h_T^{-1}\|\nu\|_{L^{q'}(D_T)}\nonumber \\ &\leqslant c\sum_{T\in\cT_{k}}\eta_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)h_T^{-1}\|\mu-I_k\mu\|_{L^{q'}(D_T)}\nonumber \\ &= c\big(\sum_{T\in\cT_{k}\setminus\cT_l^+}\eta_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)h_T\|\mu\|_{W^{2,q'}(D_T)}+\sum_{T\in\cT_l^+}\eta_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)h_T\|\mu\|_{W^{2,q'}(D_T)}\big)\nonumber \\ &\leqslant c\|h_{l}\chi_{\Omega_l^0}\|_{L^\infty(\Omega)}\big(\sum_{T\in\cT_k\setminus\cT_{l}^+}\eta^q_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\big)^{1/q}\|\mu\|_{W^{2,q'}(\Om)} +c\big(\sum_{T\in\cT_{l}^+}\eta^q_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\big)^{1/q}\|\mu\|_{W^{2,q'}(\Om)}. \nonumber \end{align*}$$

Since $\newcommand{\cT}{{\mathcal T}} \newcommand{\e}{{\rm e}} \newcommand{\bi}{\boldsymbol} \big(\sum\nolimits_{T\in\cT_k\setminus\cT_{l}^+}\eta^q_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\big){}^{1/q}\leqslant c$, see (5.16), there holds

$$\begin{align*} \newcommand{\bi}{\boldsymbol} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \newcommand{\cT}{{\mathcal T}} \displaystyle |{\rm II}|\leqslant c\|h_{l}\chi_{\Omega_l^0}\|_{L^\infty(\Omega)}\|\mu\|_{W^{2,q'}(\Om)}+c\big(\sum_{T\in\cT_{l}^+}\eta^q_{k,3}(\sigma_k^\ast,u_k^\ast,p_k^\ast,T)\big)^{1/q}\|\mu\|_{W^{2,q'}(\Om)}. \nonumber \end{align*}$$

Now by repeating the argument for the term $\rm I$, we obtain

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \label{eqn:gat_mc_aux2} {\rm II} \rightarrow 0\quad \forall\mu\in\widetilde{\mathcal{A}}\cap C^\infty(\overline{\Om}). \nonumber \end{align} \tag{ 5.18 }$$

Step iv. Take limit in preliminary variational inequality. Using (5.15) and the $\newcommand{\Om}{\Omega} H^1(\Omega)$-convergence of $\{\sigma_{k}^\ast\}_{k\geqslant0}$ in Theorem 5.2, we have for each $\newcommand{\Om}{\Omega} \mu\in\widetilde{\mathcal{A}}\cap C^\infty(\overline{\Om})$

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:gat_mc_03} \alpha'(\nabla\sigma^\ast_k,\nabla(\mu-\sigma^\ast_k)) +\alpha''(W'(\sigma_{k}^\ast),\mu-\sigma^\ast_k) \rightarrow \alpha'(\nabla\sigma^\ast_\infty,\nabla(\mu-\sigma^\ast_\infty)) +\alpha''(W'(\sigma_{\infty}^\ast),\mu-\sigma^\ast_\infty). \nonumber \end{align} \tag{ 5.19 }$$

Further, the uniform boundedness on $\{u_k^\ast\}_{k\geqslant0}$ in $\newcommand{\Om}{\Omega} H^1(\Omega)$ and the convergence of $\{\,p^\ast_{k}\}_{k\geqslant0}$ to $p_\infty^\ast$ in $\newcommand{\Om}{\Omega} H^1(\Omega)$ in theorem 5.3 yield

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle |(\mu\nabla u_k^\ast,\nabla (\,p_{k}^\ast-p_\infty^\ast))|\leqslant c\|\nabla(\,p_{k}^\ast-p_\infty^\ast)\|_{L^2(\Omega)}\rightarrow 0. \nonumber \end{align*}$$

This and theorem 5.2 imply

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \label{eqn:gat_mc_04} (\mu\nabla u_k^\ast,\nabla p_{k}^\ast)=(\mu \nabla u_k^\ast,\nabla (\,p_{k}^\ast-p_\infty^\ast))+ (\mu\nabla u_k^\ast,\nabla p_\infty^\ast) \rightarrow (\mu\nabla u^\ast_\infty,\nabla p_\infty^\ast)\quad \forall\mu\in\widetilde{\mathcal{A}}\cap C^\infty(\overline{\Om}). \nonumber \end{align} \tag{ 5.20 }$$

In the splitting

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle (\sigma_k^\ast\nabla u_{k}^\ast, \nabla p_k^\ast)-(\sigma_\infty^\ast\nabla u_{\infty}^\ast, \nabla p_\infty^\ast)&= (\sigma_k^\ast\nabla u_{k}^\ast, \nabla(\,p_k^\ast-p_\infty^\ast)) +((\sigma_k^\ast-\sigma_{\infty}^\ast)\nabla u_{k}^\ast, \nabla p_\infty^\ast)\nonumber \\ &\quad +(\sigma_{\infty}^\ast\nabla(u_{k}^\ast-u_\infty^\ast),\nabla p_\infty^\ast), \nonumber \end{align*}$$

the arguments for (5.20) directly yield

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle |(\sigma_k^\ast\nabla u_{k}^\ast, \nabla(\,p_k^\ast-p_\infty^\ast))|\rightarrow 0\quad {\rm and}\quad |(\sigma_{\infty}^\ast\nabla(u_{k}^\ast-u_\infty^\ast),\nabla p_\infty^\ast)|\rightarrow 0. \nonumber \end{align*}$$

The boundedness on $\{u_k^\ast\}_{k\geqslant0}$ in $\newcommand{\Om}{\Omega} H^1(\Omega)$, pointwise convergence of $\{\sigma_{k}^\ast\}_{k\geqslant0}$ of theorem 5.2 and Lebesgue's dominated convergence theorem imply

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle |((\sigma_k^\ast-\sigma_{\infty}^\ast)\nabla u_{k}^\ast, \nabla p_\infty^\ast)|\leqslant c\|(\sigma_k^\ast-\sigma_{\infty}^\ast)\nabla p_\infty^\ast \|_{L^2(\Omega)}\rightarrow 0. \nonumber \end{align*}$$

Hence, there holds

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{eqn:gat_mc_05} (\sigma_k^\ast\nabla u_{k}^\ast, \nabla p_k^\ast) \rightarrow (\sigma_\infty^\ast\nabla u_{\infty}^\ast, \nabla p_\infty^\ast). \nonumber \end{align} \tag{ 5.21 }$$

Now by passing both sides of (5.14) to the limit and combining the estimates (5.17)–(5.21), we obtain

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \alpha'(\nabla\sigma^\ast_\infty,\nabla(\mu-\sigma^\ast_\infty)) + \alpha''(W'(\sigma^\ast_\infty),\mu-\sigma^\ast_\infty) -(\nabla u^\ast_\infty,\nabla p_\infty^\ast(\mu-\sigma_\infty^\ast))_{L^2(\Omega)}\geqslant 0\quad \forall\mu\in\widetilde{\mathcal{A}}\cap C^\infty(\overline{\Om}). \nonumber \end{align*}$$

Step v. Density argument. By the density of $\newcommand{\Om}{\Omega} C^\infty(\overline{\Omega})$ in $\newcommand{\Om}{\Omega} H^1(\Omega)$ and the construction via a standard mollifier [19], for any $\mu\in\widetilde{\mathcal{A}}$ there exists a sequence $\newcommand{\Om}{\Omega} \{\mu_n\} \subset\widetilde{\mathcal{A}}\cap C^\infty(\overline{\Om})$ such that $\newcommand{\Om}{\Omega} \|\mu_n-\mu\|_{H^1(\Omega)}\rightarrow 0$ as $n\rightarrow\infty$. Thus, $(\nabla\sigma^\ast_\infty,\nabla \mu_n)\rightarrow(\nabla\sigma^\ast_\infty,\nabla\mu)$, $(W'(\sigma_\infty^\ast),\mu_n)\rightarrow (W'(\sigma_\infty^\ast),\mu)$, and $(\mu_n\nabla u^\ast_\infty,\nabla p_\infty^\ast)\rightarrow(\mu\nabla u^\ast_\infty,\nabla p_\infty^\ast)$, after possibly passing to a subsequence. The desired result follows from the preceding two estimates. □

Remark 5.2. The computable quantity $\newcommand{\e}{{\rm e}} \eta_{k,3}(\sigma_k^*,u_k^*,p_k^*,T)$ emerges naturally from the proof, i.e. the upper bounds on ${\rm I}$ and ${\rm II}$, which motivates its use as the a posteriori error estimator in algorithm 1.

The authors are grateful to an anonymous referee and the boarder member for the constructive comments, which have significantly improved the presentation of the paper. The work of Y Xu was partially supported by the National Natural Science Foundation of China (Grant No. 11201307), the Ministry of the Education of China through the Specialized Research Fund for the Doctoral Program of Higher Education (Grant No. 20123127120001) and the Natural Science Foundation of Shanghai (Grant No. 17ZR1420800).

Now we describe an iterative method for minimizing the energy functional

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \frac{\tilde\alpha \varepsilon}{2}\|\nabla \sigma\|_{L^2(\Omega)}^2 + \frac{\tilde\alpha}{2\varepsilon}\int_\Omega W(\sigma){\rm d}x + \frac12\|U(\sigma)-U^\delta\|^2. \nonumber \end{align*}$$

Let $p(z)=(z-c_0)(z-c_1)$. Then one linearized approximation $p_L(z,z_k)$ reads (with $\delta z=z-z_k$)

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle p_L(z,z_k) &= p(z_k) + p'(z_k)(z-z_k)\nonumber \\ &= (z_k^2-(c_0+c_1)z_k+c_0c_1) + (2z_k-c_0-c_1)\delta z. \nonumber \end{align*}$$

Upon substituting the approximation $p_L(z,z_k)$ for $p(z)$ and linearizing the forward map $U(\sigma)$, we obtain the following surrogate energy functional (with $\delta\sigma=\sigma-\sigma_k$ being the increment and $\delta U= U^\delta-U(\sigma_k)$)

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \label{eqn:surrogate} \frac{\tilde\alpha \varepsilon}{2}\|\nabla(\sigma_k+\delta\sigma)\|_{L^2(\Omega)}^2 + \frac{\tilde\alpha}{2\varepsilon}\|p(\sigma_k)+p'(\sigma_k)\delta\sigma\|_{L^2(\Omega)}^2 + \frac12\|U'(\sigma_k)\delta\sigma-\delta U\|^2. \nonumber \end{align} \tag{ A.1 }$$

The treatment of the double-well potential term $\newcommand{\Om}{\Omega} \int_\Omega W(\sigma){\rm d}x$ is in the spirit of the classical majorization–minimization algorithm in the following sense (see [56] for a detailed derivation)

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle \int_\Omega W(\sigma_k) {\rm d}x &= \int_\Omega p_L(\sigma_k,\sigma_k)^2{\rm d}x, \quad \nabla \int_\Omega W(\sigma_k) {\rm d}x = \nabla \int_\Omega p_L(\sigma_k,\sigma_k)^2{\rm d}x, \nonumber \\ \quad {\rm and} &\quad \nabla^2 \int_\Omega W(\sigma_k) {\rm d}x \leqslant \nabla^2 \int_\Omega p_L(\sigma_k,\sigma_k)^2{\rm d}x. \nonumber \end{align*}$$

This algorithm is known to have excellent numerical stability. Upon ignoring the box constraint on the conductivity $\sigma$, problem (A.1) is to find $\newcommand{\Om}{\Omega} \delta\sigma \in H^1(\Omega)$ such that

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\Om}{\Omega} \displaystyle (U'(\sigma_k)^*U'(\sigma_k)\delta\sigma,\phi) + \tilde\alpha \varepsilon& (\nabla \delta\sigma,\nabla \phi) + \frac{\tilde\alpha}{\varepsilon} (\,p'(\sigma_k)^2\delta\sigma,\phi) \nonumber \\ & = (U'(\sigma_k)^*\delta U, \phi)-\frac{\tilde\alpha}{\varepsilon}(\,p(\sigma_k)\,p'(\sigma_k),\phi)-\tilde\alpha\varepsilon(\nabla\sigma_k,\nabla\phi),\quad \forall \phi\in H^1(\Omega). \nonumber \end{align*}$$

This equation can be solved by an iterative method for the update $\delta\sigma$ (with the box constraint treated by a projection step). Note that $U'(\sigma_k)$ and $U'(\sigma_k){}^*$ can be implemented in a matrix-free manner using the standard adjoint technique. In our experiment, we employ the conjugate gradient method to solve the resulting linear systems, preconditioned by the sparse matrix corresponding to $\tilde\alpha \varepsilon (\nabla \delta\sigma,\nabla \phi) + \frac{\tilde\alpha}{\varepsilon}(\,p'(\sigma_k){}^2\delta\sigma,\phi)$.

The proof follows that in [13, 25]. By Hölder inequality and $h_T^d \leqslant |\omega_x|$ for each node $x\in T$:

$$\begin{align*} \newcommand{\e}{{\rm e}} \newcommand{\dx}{\,{\rm d}x} \displaystyle \left|\frac{1}{|\omega_x|}\int_{\omega_x}v\dx\right| \leqslant |\omega_x|^{-1/r}\|v\|_{L^r(\omega_x)}\leqslant h_T^{-d/r}\|v\|_{L^r(\omega_x)}. \nonumber \end{align*}$$

The desired L^r-stability follows from the estimate $\|\phi_x\|_{L^r(T)}\leqslant c h_T^{d/r}$, by the local quasi-uniformity of the mesh. In view of the definition (5.13), $\Pi_k \zeta = \zeta$ for any $\zeta\in\mathbb{R}$. By local inverse estimate, the L^r-stability of $\Pi_k$, standard interpolation error estimate [16] and local quasi-uniformity,

$$\begin{align} \newcommand{\e}{{\rm e}} \newcommand{\dx}{\,{\rm d}x} \displaystyle \|\nabla \Pi_k v\|_{L^r(T)}&=\inf_{\zeta\in \mathbb{R}} \|\nabla \Pi_k (v - \zeta)\|_{L^r(T)} \leqslant c h_T^{-1}\inf_{\zeta\in \mathbb{R}} \| \Pi_k (v - \zeta) \|_{L^r(T)} \nonumber \\ &\leqslant c h_T^{-1}\inf_{\zeta\in \mathbb{R}} \| v - \zeta \|_{L^r(D_T)} \leqslant c h_T^{-1} \| v - \frac{1}{|D_T|}\int_{D_T} v \dx \|_{L^r(D_T)} \leqslant c \|\nabla v \|_{L^r(D_T)}.\label{pf:int_err1} \nonumber \end{align} \tag{ B.1 }$$

Similarly,

$$\begin{align} \newcommand{\e}{{\rm e}} \displaystyle \label{pf:int_err2} \|v - \Pi_k v\|_{L^r(T)}&= \|v - \zeta - \Pi_k (v - \zeta)\|_{L^r(T)} \nonumber \\ &\leqslant c \inf_{\zeta \in \mathbb{R}} \| v - \zeta\|_{L^r(D_T)} \leqslant c h_T \|\nabla v \|_{L^r(D_T)}. \nonumber \end{align} \tag{ B.2 }$$

By the scaled trace theorem, for any $F\subset \partial T$, there holds

$$\begin{align*} \newcommand{\e}{{\rm e}} \displaystyle \|v - \Pi_k v \|_{L^r(F)} \leqslant c (h_F^{-1/r} \| v - \Pi_k v\|_{L^r(T)} + h_F^{1-1/r} \|\nabla (v - \Pi_k v) \|_{L^r(T)}). \nonumber \end{align*}$$

Then (B.1) and (B.2) complete the proof of the lemma.